Actual source code: spectraladjointassimilation.c
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: */
10: /* ------------------------------------------------------------------------
12: This program uses the one-dimensional advection-diffusion equation),
13: u_t = mu*u_xx - a 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: ------------------------------------------------------------------------- */
23: /*
24: Include "petscts.h" so that we can use TS solvers. Note that this file
25: automatically includes:
26: petscsys.h - base PETSc routines petscvec.h - vectors
27: petscmat.h - matrices
28: petscis.h - index sets petscksp.h - Krylov subspace methods
29: petscviewer.h - viewers petscpc.h - preconditioners
30: petscksp.h - linear solvers petscsnes.h - nonlinear solvers
31: */
33: #include <petsctao.h>
34: #include <petscts.h>
35: #include <petscdt.h>
36: #include <petscdraw.h>
37: #include <petscdmda.h>
39: /*
40: User-defined application context - contains data needed by the
41: application-provided call-back routines.
42: */
44: typedef struct {
45: PetscInt n; /* number of nodes */
46: PetscReal *nodes; /* GLL nodes */
47: PetscReal *weights; /* GLL weights */
48: } PetscGLL;
50: typedef struct {
51: PetscInt N; /* grid points per elements*/
52: PetscInt E; /* number of elements */
53: PetscReal tol_L2,tol_max; /* error norms */
54: PetscInt steps; /* number of timesteps */
55: PetscReal Tend; /* endtime */
56: PetscReal mu; /* viscosity */
57: PetscReal a; /* advection speed */
58: PetscReal L; /* total length of domain */
59: PetscReal Le;
60: PetscReal Tadj;
61: } PetscParam;
63: typedef struct {
64: Vec reference; /* desired end state */
65: Vec grid; /* total grid */
66: Vec grad;
67: Vec ic;
68: Vec curr_sol;
69: Vec joe;
70: Vec true_solution; /* actual initial conditions for the final solution */
71: } PetscData;
73: typedef struct {
74: Vec grid; /* total grid */
75: Vec mass; /* mass matrix for total integration */
76: Mat stiff; /* stifness matrix */
77: Mat advec;
78: Mat keptstiff;
79: PetscGLL gll;
80: } PetscSEMOperators;
82: typedef struct {
83: DM da; /* distributed array data structure */
84: PetscSEMOperators SEMop;
85: PetscParam param;
86: PetscData dat;
87: TS ts;
88: PetscReal initial_dt;
89: PetscReal *solutioncoefficients;
90: PetscInt ncoeff;
91: } AppCtx;
93: /*
94: User-defined routines
95: */
96: extern PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
97: extern PetscErrorCode RHSLaplacian(TS,PetscReal,Vec,Mat,Mat,void*);
98: extern PetscErrorCode RHSAdvection(TS,PetscReal,Vec,Mat,Mat,void*);
99: extern PetscErrorCode InitialConditions(Vec,AppCtx*);
100: extern PetscErrorCode ComputeReference(TS,PetscReal,Vec,AppCtx*);
101: extern PetscErrorCode MonitorError(Tao,void*);
102: extern PetscErrorCode MonitorDestroy(void**);
103: extern PetscErrorCode ComputeSolutionCoefficients(AppCtx*);
104: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
105: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
107: int main(int argc,char **argv)
108: {
109: AppCtx appctx; /* user-defined application context */
110: Tao tao;
111: Vec u; /* approximate solution vector */
112: PetscInt i, xs, xm, ind, j, lenglob;
113: PetscReal x, *wrk_ptr1, *wrk_ptr2;
114: MatNullSpace nsp;
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Initialize program and set problem parameters
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: PetscInitialize(&argc,&argv,(char*)0,help);
122: /*initialize parameters */
123: appctx.param.N = 10; /* order of the spectral element */
124: appctx.param.E = 8; /* number of elements */
125: appctx.param.L = 1.0; /* length of the domain */
126: appctx.param.mu = 0.00001; /* diffusion coefficient */
127: appctx.param.a = 0.0; /* advection speed */
128: appctx.initial_dt = 1e-4;
129: appctx.param.steps = PETSC_MAX_INT;
130: appctx.param.Tend = 0.01;
131: appctx.ncoeff = 2;
133: PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL);
134: PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL);
135: PetscOptionsGetInt(NULL,NULL,"-ncoeff",&appctx.ncoeff,NULL);
136: PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL);
137: PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL);
138: PetscOptionsGetReal(NULL,NULL,"-a",&appctx.param.a,NULL);
139: appctx.param.Le = appctx.param.L/appctx.param.E;
141: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: Create GLL data structures
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: PetscMalloc2(appctx.param.N,&appctx.SEMop.gll.nodes,appctx.param.N,&appctx.SEMop.gll.weights);
145: PetscDTGaussLobattoLegendreQuadrature(appctx.param.N,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
146: appctx.SEMop.gll.n = appctx.param.N;
147: lenglob = appctx.param.E*(appctx.param.N-1);
149: /*
150: Create distributed array (DMDA) to manage parallel grid and vectors
151: and to set up the ghost point communication pattern. There are E*(Nl-1)+1
152: total grid values spread equally among all the processors, except first and last
153: */
155: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);
156: DMSetFromOptions(appctx.da);
157: DMSetUp(appctx.da);
159: /*
160: Extract global and local vectors from DMDA; we use these to store the
161: approximate solution. Then duplicate these for remaining vectors that
162: have the same types.
163: */
165: DMCreateGlobalVector(appctx.da,&u);
166: VecDuplicate(u,&appctx.dat.ic);
167: VecDuplicate(u,&appctx.dat.true_solution);
168: VecDuplicate(u,&appctx.dat.reference);
169: VecDuplicate(u,&appctx.SEMop.grid);
170: VecDuplicate(u,&appctx.SEMop.mass);
171: VecDuplicate(u,&appctx.dat.curr_sol);
172: VecDuplicate(u,&appctx.dat.joe);
174: DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
175: DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
176: DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
178: /* Compute function over the locally owned part of the grid */
180: xs=xs/(appctx.param.N-1);
181: xm=xm/(appctx.param.N-1);
183: /*
184: Build total grid and mass over entire mesh (multi-elemental)
185: */
187: for (i=xs; i<xs+xm; i++) {
188: for (j=0; j<appctx.param.N-1; j++) {
189: x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i;
190: ind=i*(appctx.param.N-1)+j;
191: wrk_ptr1[ind]=x;
192: wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
193: if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
194: }
195: }
196: DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
197: DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
199: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200: Create matrix data structure; set matrix evaluation routine.
201: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
202: DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);
203: DMCreateMatrix(appctx.da,&appctx.SEMop.stiff);
204: DMCreateMatrix(appctx.da,&appctx.SEMop.advec);
206: /*
207: For linear problems with a time-dependent f(u,t) in the equation
208: u_t = f(u,t), the user provides the discretized right-hand-side
209: as a time-dependent matrix.
210: */
211: RHSLaplacian(appctx.ts,0.0,u,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx);
212: RHSAdvection(appctx.ts,0.0,u,appctx.SEMop.advec,appctx.SEMop.advec,&appctx);
213: MatAXPY(appctx.SEMop.stiff,-1.0,appctx.SEMop.advec,DIFFERENT_NONZERO_PATTERN);
214: MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff);
216: /* attach the null space to the matrix, this probably is not needed but does no harm */
217: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
218: MatSetNullSpace(appctx.SEMop.stiff,nsp);
219: MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL);
220: MatNullSpaceDestroy(&nsp);
222: /* Create the TS solver that solves the ODE and its adjoint; set its options */
223: TSCreate(PETSC_COMM_WORLD,&appctx.ts);
224: TSSetSolutionFunction(appctx.ts,(PetscErrorCode (*)(TS,PetscReal,Vec, void *))ComputeReference,&appctx);
225: TSSetProblemType(appctx.ts,TS_LINEAR);
226: TSSetType(appctx.ts,TSRK);
227: TSSetDM(appctx.ts,appctx.da);
228: TSSetTime(appctx.ts,0.0);
229: TSSetTimeStep(appctx.ts,appctx.initial_dt);
230: TSSetMaxSteps(appctx.ts,appctx.param.steps);
231: TSSetMaxTime(appctx.ts,appctx.param.Tend);
232: TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
233: TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL);
234: TSSetFromOptions(appctx.ts);
235: /* Need to save initial timestep user may have set with -ts_dt so it can be reset for each new TSSolve() */
236: TSGetTimeStep(appctx.ts,&appctx.initial_dt);
237: TSSetRHSFunction(appctx.ts,NULL,TSComputeRHSFunctionLinear,&appctx);
238: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,TSComputeRHSJacobianConstant,&appctx);
239: /* TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
240: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx); */
242: /* Set random initial conditions as initial guess, compute analytic reference solution and analytic (true) initial conditions */
243: ComputeSolutionCoefficients(&appctx);
244: InitialConditions(appctx.dat.ic,&appctx);
245: ComputeReference(appctx.ts,appctx.param.Tend,appctx.dat.reference,&appctx);
246: ComputeReference(appctx.ts,0.0,appctx.dat.true_solution,&appctx);
248: /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */
249: TSSetSaveTrajectory(appctx.ts);
250: TSSetFromOptions(appctx.ts);
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: TaoSetSolution(tao,appctx.dat.ic);
257: /* Set routine for function and gradient evaluation */
258: TaoSetObjectiveAndGradient(tao,NULL,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: PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
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 0;
289: }
291: /*
292: Computes the coefficients for the analytic solution to the PDE
293: */
294: PetscErrorCode ComputeSolutionCoefficients(AppCtx *appctx)
295: {
296: PetscRandom rand;
297: PetscInt i;
299: PetscMalloc1(appctx->ncoeff,&appctx->solutioncoefficients);
300: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
301: PetscRandomSetInterval(rand,.9,1.0);
302: for (i=0; i<appctx->ncoeff; i++) {
303: PetscRandomGetValue(rand,&appctx->solutioncoefficients[i]);
304: }
305: PetscRandomDestroy(&rand);
306: return 0;
307: }
309: /* --------------------------------------------------------------------- */
310: /*
311: InitialConditions - Computes the (random) initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()
313: Input Parameter:
314: u - uninitialized solution vector (global)
315: appctx - user-defined application context
317: Output Parameter:
318: u - vector with solution at initial time (global)
319: */
320: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
321: {
322: PetscScalar *s;
323: const PetscScalar *xg;
324: PetscInt i,j,lenglob;
325: PetscReal sum,val;
326: PetscRandom rand;
328: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
329: PetscRandomSetInterval(rand,.9,1.0);
330: DMDAVecGetArray(appctx->da,u,&s);
331: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
332: lenglob = appctx->param.E*(appctx->param.N-1);
333: for (i=0; i<lenglob; i++) {
334: s[i]= 0;
335: for (j=0; j<appctx->ncoeff; j++) {
336: PetscRandomGetValue(rand,&val);
337: s[i] += val*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
338: }
339: }
340: PetscRandomDestroy(&rand);
341: DMDAVecRestoreArray(appctx->da,u,&s);
342: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
343: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
344: VecSum(u,&sum);
345: VecShift(u,-sum/lenglob);
346: return 0;
347: }
349: /*
350: TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function.
352: InitialConditions() computes the initial conditions for the beginning of the Tao iterations
354: Input Parameter:
355: u - uninitialized solution vector (global)
356: appctx - user-defined application context
358: Output Parameter:
359: u - vector with solution at initial time (global)
360: */
361: PetscErrorCode TrueSolution(Vec u,AppCtx *appctx)
362: {
363: PetscScalar *s;
364: const PetscScalar *xg;
365: PetscInt i,j,lenglob;
366: PetscReal sum;
368: DMDAVecGetArray(appctx->da,u,&s);
369: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
370: lenglob = appctx->param.E*(appctx->param.N-1);
371: for (i=0; i<lenglob; i++) {
372: s[i]= 0;
373: for (j=0; j<appctx->ncoeff; j++) {
374: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
375: }
376: }
377: DMDAVecRestoreArray(appctx->da,u,&s);
378: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
379: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
380: VecSum(u,&sum);
381: VecShift(u,-sum/lenglob);
382: return 0;
383: }
384: /* --------------------------------------------------------------------- */
385: /*
386: Sets the desired profile for the final end time
388: Input Parameters:
389: t - final time
390: obj - vector storing the desired profile
391: appctx - user-defined application context
393: */
394: PetscErrorCode ComputeReference(TS ts,PetscReal t,Vec obj,AppCtx *appctx)
395: {
396: PetscScalar *s,tc;
397: const PetscScalar *xg;
398: PetscInt i, j,lenglob;
400: DMDAVecGetArray(appctx->da,obj,&s);
401: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
402: lenglob = appctx->param.E*(appctx->param.N-1);
403: for (i=0; i<lenglob; i++) {
404: s[i]= 0;
405: for (j=0; j<appctx->ncoeff; j++) {
406: tc = -appctx->param.mu*(j+1)*(j+1)*4.0*PETSC_PI*PETSC_PI*t;
407: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*(xg[i] + appctx->param.a*t))*PetscExpReal(tc);
408: }
409: }
410: DMDAVecRestoreArray(appctx->da,obj,&s);
411: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
412: return 0;
413: }
415: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
416: {
417: AppCtx *appctx = (AppCtx*)ctx;
419: MatMult(appctx->SEMop.keptstiff,globalin,globalout);
420: return 0;
421: }
423: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx)
424: {
425: AppCtx *appctx = (AppCtx*)ctx;
427: MatCopy(appctx->SEMop.keptstiff,A,DIFFERENT_NONZERO_PATTERN);
428: return 0;
429: }
431: /* --------------------------------------------------------------------- */
433: /*
434: RHSLaplacian - matrix for diffusion
436: Input Parameters:
437: ts - the TS context
438: t - current time (ignored)
439: X - current solution (ignored)
440: dummy - optional user-defined context, as set by TSetRHSJacobian()
442: Output Parameters:
443: AA - Jacobian matrix
444: BB - optionally different matrix from which the preconditioner is built
445: str - flag indicating matrix structure
447: Scales by the inverse of the mass matrix (perhaps that should be pulled out)
449: */
450: PetscErrorCode RHSLaplacian(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
451: {
452: PetscReal **temp;
453: PetscReal vv;
454: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
455: PetscInt i,xs,xn,l,j;
456: PetscInt *rowsDM;
458: /*
459: Creates the element stiffness matrix for the given gll
460: */
461: PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
463: /* scale by the size of the element */
464: for (i=0; i<appctx->param.N; i++) {
465: vv=-appctx->param.mu*2.0/appctx->param.Le;
466: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
467: }
469: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
470: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
473: xs = xs/(appctx->param.N-1);
474: xn = xn/(appctx->param.N-1);
476: PetscMalloc1(appctx->param.N,&rowsDM);
477: /*
478: loop over local elements
479: */
480: for (j=xs; j<xs+xn; j++) {
481: for (l=0; l<appctx->param.N; l++) {
482: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
483: }
484: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
485: }
486: PetscFree(rowsDM);
487: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
488: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
489: VecReciprocal(appctx->SEMop.mass);
490: MatDiagonalScale(A,appctx->SEMop.mass,0);
491: VecReciprocal(appctx->SEMop.mass);
493: PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
494: return 0;
495: }
497: /*
498: Almost identical to Laplacian
500: Note that the element matrix is NOT scaled by the size of element like the Laplacian term.
501: */
502: PetscErrorCode RHSAdvection(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
503: {
504: PetscReal **temp;
505: PetscReal vv;
506: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
507: PetscInt i,xs,xn,l,j;
508: PetscInt *rowsDM;
510: /*
511: Creates the element stiffness matrix for the given gll
512: */
513: PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
515: /* scale by the size of the element */
516: for (i=0; i<appctx->param.N; i++) {
517: vv = -appctx->param.a;
518: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
519: }
521: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
522: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
525: xs = xs/(appctx->param.N-1);
526: xn = xn/(appctx->param.N-1);
528: PetscMalloc1(appctx->param.N,&rowsDM);
529: /*
530: loop over local elements
531: */
532: for (j=xs; j<xs+xn; j++) {
533: for (l=0; l<appctx->param.N; l++) {
534: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
535: }
536: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
537: }
538: PetscFree(rowsDM);
539: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
540: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
541: VecReciprocal(appctx->SEMop.mass);
542: MatDiagonalScale(A,appctx->SEMop.mass,0);
543: VecReciprocal(appctx->SEMop.mass);
545: PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
546: return 0;
547: }
549: /* ------------------------------------------------------------------ */
550: /*
551: FormFunctionGradient - Evaluates the function and corresponding gradient.
553: Input Parameters:
554: tao - the Tao context
555: ic - the input vector
556: ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradient()
558: Output Parameters:
559: f - the newly evaluated function
560: G - the newly evaluated gradient
562: Notes:
564: The forward equation is
565: M u_t = F(U)
566: which is converted to
567: u_t = M^{-1} F(u)
568: in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is
569: M^{-1} J
570: where J is the Jacobian of F. Now the adjoint equation is
571: M v_t = J^T v
572: but TSAdjoint does not solve this since it can only solve the transposed system for the
573: Jacobian the user provided. Hence TSAdjoint solves
574: w_t = J^T M^{-1} w (where w = M v)
575: since there is no way to indicate the mass matrix as a separate entity to TS. Thus one
576: must be careful in initializing the "adjoint equation" and using the result. This is
577: why
578: G = -2 M(u(T) - u_d)
579: below (instead of -2(u(T) - u_d)
581: */
582: PetscErrorCode FormFunctionGradient(Tao tao,Vec ic,PetscReal *f,Vec G,void *ctx)
583: {
584: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
585: Vec temp;
587: TSSetTime(appctx->ts,0.0);
588: TSSetStepNumber(appctx->ts,0);
589: TSSetTimeStep(appctx->ts,appctx->initial_dt);
590: VecCopy(ic,appctx->dat.curr_sol);
592: TSSolve(appctx->ts,appctx->dat.curr_sol);
593: VecCopy(appctx->dat.curr_sol,appctx->dat.joe);
595: /* Compute the difference between the current ODE solution and target ODE solution */
596: VecWAXPY(G,-1.0,appctx->dat.curr_sol,appctx->dat.reference);
598: /* Compute the objective/cost function */
599: VecDuplicate(G,&temp);
600: VecPointwiseMult(temp,G,G);
601: VecDot(temp,appctx->SEMop.mass,f);
602: VecDestroy(&temp);
604: /* Compute initial conditions for the adjoint integration. See Notes above */
605: VecScale(G, -2.0);
606: VecPointwiseMult(G,G,appctx->SEMop.mass);
607: TSSetCostGradients(appctx->ts,1,&G,NULL);
609: TSAdjointSolve(appctx->ts);
610: /* VecPointwiseDivide(G,G,appctx->SEMop.mass);*/
611: return 0;
612: }
614: PetscErrorCode MonitorError(Tao tao,void *ctx)
615: {
616: AppCtx *appctx = (AppCtx*)ctx;
617: Vec temp,grad;
618: PetscReal nrm;
619: PetscInt its;
620: PetscReal fct,gnorm;
622: VecDuplicate(appctx->dat.ic,&temp);
623: VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution);
624: VecPointwiseMult(temp,temp,temp);
625: VecDot(temp,appctx->SEMop.mass,&nrm);
626: nrm = PetscSqrtReal(nrm);
627: TaoGetGradient(tao,&grad,NULL,NULL);
628: VecPointwiseMult(temp,temp,temp);
629: VecDot(temp,appctx->SEMop.mass,&gnorm);
630: gnorm = PetscSqrtReal(gnorm);
631: VecDestroy(&temp);
632: TaoGetIterationNumber(tao,&its);
633: TaoGetSolutionStatus(tao,NULL,&fct,NULL,NULL,NULL,NULL);
634: if (!its) {
635: PetscPrintf(PETSC_COMM_WORLD,"%% Iteration Error Objective Gradient-norm\n");
636: PetscPrintf(PETSC_COMM_WORLD,"history = [\n");
637: }
638: PetscPrintf(PETSC_COMM_WORLD,"%3D %g %g %g\n",its,(double)nrm,(double)fct,(double)gnorm);
639: return 0;
640: }
642: PetscErrorCode MonitorDestroy(void **ctx)
643: {
644: PetscPrintf(PETSC_COMM_WORLD,"];\n");
645: return 0;
646: }
648: /*TEST
650: build:
651: requires: !complex
653: test:
654: requires: !single
655: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
657: test:
658: suffix: cn
659: requires: !single
660: args: -ts_type cn -ts_dt .003 -pc_type lu -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
662: test:
663: suffix: 2
664: requires: !single
665: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -a .1 -tao_bqnls_mat_lmvm_scale_type none
667: TEST*/