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: */
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 <petscdt.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; /* number of nodes */
53: PetscReal *nodes; /* GLL nodes */
54: PetscReal *weights; /* GLL weights */
55: } PetscGLL;
57: typedef struct {
58: PetscInt N; /* grid points per elements*/
59: PetscInt E; /* number of elements */
60: PetscReal tol_L2,tol_max; /* error norms */
61: PetscInt steps; /* number of timesteps */
62: PetscReal Tend; /* endtime */
63: PetscReal mu; /* viscosity */
64: PetscReal a; /* advection speed */
65: PetscReal L; /* total length of domain */
66: PetscReal Le;
67: PetscReal Tadj;
68: } PetscParam;
70: typedef struct {
71: Vec reference; /* desired end state */
72: Vec grid; /* total grid */
73: Vec grad;
74: Vec ic;
75: Vec curr_sol;
76: Vec joe;
77: Vec true_solution; /* actual initial conditions for the final solution */
78: } PetscData;
80: typedef struct {
81: Vec grid; /* total grid */
82: Vec mass; /* mass matrix for total integration */
83: Mat stiff; /* stifness matrix */
84: Mat advec;
85: Mat keptstiff;
86: PetscGLL gll;
87: } PetscSEMOperators;
89: typedef struct {
90: DM da; /* distributed array data structure */
91: PetscSEMOperators SEMop;
92: PetscParam param;
93: PetscData dat;
94: TS ts;
95: PetscReal initial_dt;
96: PetscReal *solutioncoefficients;
97: PetscInt ncoeff;
98: } AppCtx;
100: /*
101: User-defined routines
102: */
103: extern PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
104: extern PetscErrorCode RHSLaplacian(TS,PetscReal,Vec,Mat,Mat,void*);
105: extern PetscErrorCode RHSAdvection(TS,PetscReal,Vec,Mat,Mat,void*);
106: extern PetscErrorCode InitialConditions(Vec,AppCtx*);
107: extern PetscErrorCode ComputeReference(TS,PetscReal,Vec,AppCtx*);
108: extern PetscErrorCode MonitorError(Tao,void*);
109: extern PetscErrorCode MonitorDestroy(void**);
110: extern PetscErrorCode ComputeSolutionCoefficients(AppCtx*);
111: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
112: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
114: int main(int argc,char **argv)
115: {
116: AppCtx appctx; /* user-defined application context */
117: Tao tao;
118: Vec u; /* approximate solution vector */
120: PetscInt i, xs, xm, ind, j, lenglob;
121: PetscReal x, *wrk_ptr1, *wrk_ptr2;
122: MatNullSpace nsp;
124: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: Initialize program and set problem parameters
126: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
131: /*initialize parameters */
132: appctx.param.N = 10; /* order of the spectral element */
133: appctx.param.E = 8; /* number of elements */
134: appctx.param.L = 1.0; /* length of the domain */
135: appctx.param.mu = 0.00001; /* diffusion coefficient */
136: appctx.param.a = 0.0; /* advection speed */
137: appctx.initial_dt = 1e-4;
138: appctx.param.steps = PETSC_MAX_INT;
139: appctx.param.Tend = 0.01;
140: appctx.ncoeff = 2;
142: PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL);
143: PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL);
144: PetscOptionsGetInt(NULL,NULL,"-ncoeff",&appctx.ncoeff,NULL);
145: PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL);
146: PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL);
147: PetscOptionsGetReal(NULL,NULL,"-a",&appctx.param.a,NULL);
148: appctx.param.Le = appctx.param.L/appctx.param.E;
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Create GLL data structures
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: PetscMalloc2(appctx.param.N,&appctx.SEMop.gll.nodes,appctx.param.N,&appctx.SEMop.gll.weights);
154: PetscDTGaussLobattoLegendreQuadrature(appctx.param.N,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
155: appctx.SEMop.gll.n = appctx.param.N;
156: lenglob = appctx.param.E*(appctx.param.N-1);
158: /*
159: Create distributed array (DMDA) to manage parallel grid and vectors
160: and to set up the ghost point communication pattern. There are E*(Nl-1)+1
161: total grid values spread equally among all the processors, except first and last
162: */
164: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);
165: DMSetFromOptions(appctx.da);
166: DMSetUp(appctx.da);
168: /*
169: Extract global and local vectors from DMDA; we use these to store the
170: approximate solution. Then duplicate these for remaining vectors that
171: have the same types.
172: */
174: DMCreateGlobalVector(appctx.da,&u);
175: VecDuplicate(u,&appctx.dat.ic);
176: VecDuplicate(u,&appctx.dat.true_solution);
177: VecDuplicate(u,&appctx.dat.reference);
178: VecDuplicate(u,&appctx.SEMop.grid);
179: VecDuplicate(u,&appctx.SEMop.mass);
180: VecDuplicate(u,&appctx.dat.curr_sol);
181: VecDuplicate(u,&appctx.dat.joe);
183: DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
184: DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
185: DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
187: /* Compute function over the locally owned part of the grid */
189: xs=xs/(appctx.param.N-1);
190: xm=xm/(appctx.param.N-1);
192: /*
193: Build total grid and mass over entire mesh (multi-elemental)
194: */
196: for (i=xs; i<xs+xm; i++) {
197: for (j=0; j<appctx.param.N-1; j++) {
198: x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i;
199: ind=i*(appctx.param.N-1)+j;
200: wrk_ptr1[ind]=x;
201: wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
202: if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
203: }
204: }
205: DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
206: DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
208: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
209: Create matrix data structure; set matrix evaluation routine.
210: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
211: DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);
212: DMCreateMatrix(appctx.da,&appctx.SEMop.stiff);
213: DMCreateMatrix(appctx.da,&appctx.SEMop.advec);
215: /*
216: For linear problems with a time-dependent f(u,t) in the equation
217: u_t = f(u,t), the user provides the discretized right-hand-side
218: as a time-dependent matrix.
219: */
220: RHSLaplacian(appctx.ts,0.0,u,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx);
221: RHSAdvection(appctx.ts,0.0,u,appctx.SEMop.advec,appctx.SEMop.advec,&appctx);
222: MatAXPY(appctx.SEMop.stiff,-1.0,appctx.SEMop.advec,DIFFERENT_NONZERO_PATTERN);
223: MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff);
225: /* attach the null space to the matrix, this probably is not needed but does no harm */
226: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
227: MatSetNullSpace(appctx.SEMop.stiff,nsp);
228: MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL);
229: MatNullSpaceDestroy(&nsp);
231: /* Create the TS solver that solves the ODE and its adjoint; set its options */
232: TSCreate(PETSC_COMM_WORLD,&appctx.ts);
233: TSSetSolutionFunction(appctx.ts,(PetscErrorCode (*)(TS,PetscReal,Vec, void *))ComputeReference,&appctx);
234: TSSetProblemType(appctx.ts,TS_LINEAR);
235: TSSetType(appctx.ts,TSRK);
236: TSSetDM(appctx.ts,appctx.da);
237: TSSetTime(appctx.ts,0.0);
238: TSSetTimeStep(appctx.ts,appctx.initial_dt);
239: TSSetMaxSteps(appctx.ts,appctx.param.steps);
240: TSSetMaxTime(appctx.ts,appctx.param.Tend);
241: TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
242: TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL);
243: TSSetFromOptions(appctx.ts);
244: /* Need to save initial timestep user may have set with -ts_dt so it can be reset for each new TSSolve() */
245: TSGetTimeStep(appctx.ts,&appctx.initial_dt);
246: TSSetRHSFunction(appctx.ts,NULL,TSComputeRHSFunctionLinear,&appctx);
247: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,TSComputeRHSJacobianConstant,&appctx);
248: /* TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
249: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx); */
251: /* Set random initial conditions as initial guess, compute analytic reference solution and analytic (true) initial conditions */
252: ComputeSolutionCoefficients(&appctx);
253: InitialConditions(appctx.dat.ic,&appctx);
254: ComputeReference(appctx.ts,appctx.param.Tend,appctx.dat.reference,&appctx);
255: ComputeReference(appctx.ts,0.0,appctx.dat.true_solution,&appctx);
257: /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */
258: TSSetSaveTrajectory(appctx.ts);
259: TSSetFromOptions(appctx.ts);
261: /* Create TAO solver and set desired solution method */
262: TaoCreate(PETSC_COMM_WORLD,&tao);
263: TaoSetMonitor(tao,MonitorError,&appctx,MonitorDestroy);
264: TaoSetType(tao,TAOBQNLS);
265: TaoSetInitialVector(tao,appctx.dat.ic);
266: /* Set routine for function and gradient evaluation */
267: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&appctx);
268: /* Check for any TAO command line options */
269: TaoSetTolerances(tao,1e-8,PETSC_DEFAULT,PETSC_DEFAULT);
270: TaoSetFromOptions(tao);
271: TaoSolve(tao);
273: TaoDestroy(&tao);
274: PetscFree(appctx.solutioncoefficients);
275: MatDestroy(&appctx.SEMop.advec);
276: MatDestroy(&appctx.SEMop.stiff);
277: MatDestroy(&appctx.SEMop.keptstiff);
278: VecDestroy(&u);
279: VecDestroy(&appctx.dat.ic);
280: VecDestroy(&appctx.dat.joe);
281: VecDestroy(&appctx.dat.true_solution);
282: VecDestroy(&appctx.dat.reference);
283: VecDestroy(&appctx.SEMop.grid);
284: VecDestroy(&appctx.SEMop.mass);
285: VecDestroy(&appctx.dat.curr_sol);
286: PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
287: DMDestroy(&appctx.da);
288: TSDestroy(&appctx.ts);
290: /*
291: Always call PetscFinalize() before exiting a program. This routine
292: - finalizes the PETSc libraries as well as MPI
293: - provides summary and diagnostic information if certain runtime
294: options are chosen (e.g., -log_summary).
295: */
296: PetscFinalize();
297: return ierr;
298: }
300: /*
301: Computes the coefficients for the analytic solution to the PDE
302: */
303: PetscErrorCode ComputeSolutionCoefficients(AppCtx *appctx)
304: {
305: PetscErrorCode ierr;
306: PetscRandom rand;
307: PetscInt i;
310: PetscMalloc1(appctx->ncoeff,&appctx->solutioncoefficients);
311: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
312: PetscRandomSetInterval(rand,.9,1.0);
313: for (i=0; i<appctx->ncoeff; i++) {
314: PetscRandomGetValue(rand,&appctx->solutioncoefficients[i]);
315: }
316: PetscRandomDestroy(&rand);
317: return(0);
318: }
320: /* --------------------------------------------------------------------- */
321: /*
322: InitialConditions - Computes the (random) initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()
324: Input Parameter:
325: u - uninitialized solution vector (global)
326: appctx - user-defined application context
328: Output Parameter:
329: u - vector with solution at initial time (global)
330: */
331: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
332: {
333: PetscScalar *s;
334: const PetscScalar *xg;
335: PetscErrorCode ierr;
336: PetscInt i,j,lenglob;
337: PetscReal sum,val;
338: PetscRandom rand;
341: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
342: PetscRandomSetInterval(rand,.9,1.0);
343: DMDAVecGetArray(appctx->da,u,&s);
344: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
345: lenglob = appctx->param.E*(appctx->param.N-1);
346: for (i=0; i<lenglob; i++) {
347: s[i]= 0;
348: for (j=0; j<appctx->ncoeff; j++) {
349: PetscRandomGetValue(rand,&val);
350: s[i] += val*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
351: }
352: }
353: PetscRandomDestroy(&rand);
354: DMDAVecRestoreArray(appctx->da,u,&s);
355: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
356: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
357: VecSum(u,&sum);
358: VecShift(u,-sum/lenglob);
359: return(0);
360: }
362: /*
363: TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function.
365: InitialConditions() computes the initial conditions for the beginning of the Tao iterations
367: Input Parameter:
368: u - uninitialized solution vector (global)
369: appctx - user-defined application context
371: Output Parameter:
372: u - vector with solution at initial time (global)
373: */
374: PetscErrorCode TrueSolution(Vec u,AppCtx *appctx)
375: {
376: PetscScalar *s;
377: const PetscScalar *xg;
378: PetscErrorCode ierr;
379: PetscInt i,j,lenglob;
380: PetscReal sum;
383: DMDAVecGetArray(appctx->da,u,&s);
384: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
385: lenglob = appctx->param.E*(appctx->param.N-1);
386: for (i=0; i<lenglob; i++) {
387: s[i]= 0;
388: for (j=0; j<appctx->ncoeff; j++) {
389: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
390: }
391: }
392: DMDAVecRestoreArray(appctx->da,u,&s);
393: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
394: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
395: VecSum(u,&sum);
396: VecShift(u,-sum/lenglob);
397: return(0);
398: }
399: /* --------------------------------------------------------------------- */
400: /*
401: Sets the desired profile for the final end time
403: Input Parameters:
404: t - final time
405: obj - vector storing the desired profile
406: appctx - user-defined application context
408: */
409: PetscErrorCode ComputeReference(TS ts,PetscReal t,Vec obj,AppCtx *appctx)
410: {
411: PetscScalar *s,tc;
412: const PetscScalar *xg;
413: PetscErrorCode ierr;
414: PetscInt i, j,lenglob;
417: DMDAVecGetArray(appctx->da,obj,&s);
418: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
419: lenglob = appctx->param.E*(appctx->param.N-1);
420: for (i=0; i<lenglob; i++) {
421: s[i]= 0;
422: for (j=0; j<appctx->ncoeff; j++) {
423: tc = -appctx->param.mu*(j+1)*(j+1)*4.0*PETSC_PI*PETSC_PI*t;
424: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*(xg[i] + appctx->param.a*t))*PetscExpReal(tc);
425: }
426: }
427: DMDAVecRestoreArray(appctx->da,obj,&s);
428: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
429: return(0);
430: }
432: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
433: {
435: AppCtx *appctx = (AppCtx*)ctx;
438: MatMult(appctx->SEMop.keptstiff,globalin,globalout);
439: return(0);
440: }
442: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx)
443: {
445: AppCtx *appctx = (AppCtx*)ctx;
448: MatCopy(appctx->SEMop.keptstiff,A,DIFFERENT_NONZERO_PATTERN);
449: return(0);
450: }
452: /* --------------------------------------------------------------------- */
454: /*
455: RHSLaplacian - matrix for diffusion
457: Input Parameters:
458: ts - the TS context
459: t - current time (ignored)
460: X - current solution (ignored)
461: dummy - optional user-defined context, as set by TSetRHSJacobian()
463: Output Parameters:
464: AA - Jacobian matrix
465: BB - optionally different matrix from which the preconditioner is built
466: str - flag indicating matrix structure
468: Scales by the inverse of the mass matrix (perhaps that should be pulled out)
470: */
471: PetscErrorCode RHSLaplacian(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
472: {
473: PetscReal **temp;
474: PetscReal vv;
475: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
477: PetscInt i,xs,xn,l,j;
478: PetscInt *rowsDM;
481: /*
482: Creates the element stiffness matrix for the given gll
483: */
484: PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
486: /* scale by the size of the element */
487: for (i=0; i<appctx->param.N; i++) {
488: vv=-appctx->param.mu*2.0/appctx->param.Le;
489: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
490: }
492: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
493: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
495: if (appctx->param.N-1 < 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Polynomial order must be at least 2");
496: xs = xs/(appctx->param.N-1);
497: xn = xn/(appctx->param.N-1);
499: PetscMalloc1(appctx->param.N,&rowsDM);
500: /*
501: loop over local elements
502: */
503: for (j=xs; j<xs+xn; j++) {
504: for (l=0; l<appctx->param.N; l++) {
505: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
506: }
507: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
508: }
509: PetscFree(rowsDM);
510: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
511: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
512: VecReciprocal(appctx->SEMop.mass);
513: MatDiagonalScale(A,appctx->SEMop.mass,0);
514: VecReciprocal(appctx->SEMop.mass);
516: PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
517: return(0);
518: }
520: /*
521: Almost identical to Laplacian
523: Note that the element matrix is NOT scaled by the size of element like the Laplacian term.
524: */
525: PetscErrorCode RHSAdvection(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
526: {
527: PetscReal **temp;
528: PetscReal vv;
529: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
531: PetscInt i,xs,xn,l,j;
532: PetscInt *rowsDM;
535: /*
536: Creates the element stiffness matrix for the given gll
537: */
538: PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
540: /* scale by the size of the element */
541: for (i=0; i<appctx->param.N; i++) {
542: vv = -appctx->param.a;
543: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
544: }
546: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
547: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
549: if (appctx->param.N-1 < 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Polynomial order must be at least 2");
550: xs = xs/(appctx->param.N-1);
551: xn = xn/(appctx->param.N-1);
553: PetscMalloc1(appctx->param.N,&rowsDM);
554: /*
555: loop over local elements
556: */
557: for (j=xs; j<xs+xn; j++) {
558: for (l=0; l<appctx->param.N; l++) {
559: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
560: }
561: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
562: }
563: PetscFree(rowsDM);
564: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
565: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
566: VecReciprocal(appctx->SEMop.mass);
567: MatDiagonalScale(A,appctx->SEMop.mass,0);
568: VecReciprocal(appctx->SEMop.mass);
570: PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
571: return(0);
572: }
574: /* ------------------------------------------------------------------ */
575: /*
576: FormFunctionGradient - Evaluates the function and corresponding gradient.
578: Input Parameters:
579: tao - the Tao context
580: ic - the input vector
581: ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradientRoutine()
583: Output Parameters:
584: f - the newly evaluated function
585: G - the newly evaluated gradient
587: Notes:
589: The forward equation is
590: M u_t = F(U)
591: which is converted to
592: u_t = M^{-1} F(u)
593: in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is
594: M^{-1} J
595: where J is the Jacobian of F. Now the adjoint equation is
596: M v_t = J^T v
597: but TSAdjoint does not solve this since it can only solve the transposed system for the
598: Jacobian the user provided. Hence TSAdjoint solves
599: w_t = J^T M^{-1} w (where w = M v)
600: since there is no way to indicate the mass matrix as a separate entity to TS. Thus one
601: must be careful in initializing the "adjoint equation" and using the result. This is
602: why
603: G = -2 M(u(T) - u_d)
604: below (instead of -2(u(T) - u_d)
606: */
607: PetscErrorCode FormFunctionGradient(Tao tao,Vec ic,PetscReal *f,Vec G,void *ctx)
608: {
609: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
610: PetscErrorCode ierr;
611: Vec temp;
614: TSSetTime(appctx->ts,0.0);
615: TSSetStepNumber(appctx->ts,0);
616: TSSetTimeStep(appctx->ts,appctx->initial_dt);
617: VecCopy(ic,appctx->dat.curr_sol);
619: TSSolve(appctx->ts,appctx->dat.curr_sol);
620: VecCopy(appctx->dat.curr_sol,appctx->dat.joe);
622: /* Compute the difference between the current ODE solution and target ODE solution */
623: VecWAXPY(G,-1.0,appctx->dat.curr_sol,appctx->dat.reference);
625: /* Compute the objective/cost function */
626: VecDuplicate(G,&temp);
627: VecPointwiseMult(temp,G,G);
628: VecDot(temp,appctx->SEMop.mass,f);
629: VecDestroy(&temp);
631: /* Compute initial conditions for the adjoint integration. See Notes above */
632: VecScale(G, -2.0);
633: VecPointwiseMult(G,G,appctx->SEMop.mass);
634: TSSetCostGradients(appctx->ts,1,&G,NULL);
636: TSAdjointSolve(appctx->ts);
637: /* VecPointwiseDivide(G,G,appctx->SEMop.mass);*/
638: return(0);
639: }
641: PetscErrorCode MonitorError(Tao tao,void *ctx)
642: {
643: AppCtx *appctx = (AppCtx*)ctx;
644: Vec temp,grad;
645: PetscReal nrm;
647: PetscInt its;
648: PetscReal fct,gnorm;
651: VecDuplicate(appctx->dat.ic,&temp);
652: VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution);
653: VecPointwiseMult(temp,temp,temp);
654: VecDot(temp,appctx->SEMop.mass,&nrm);
655: nrm = PetscSqrtReal(nrm);
656: TaoGetGradientVector(tao,&grad);
657: VecPointwiseMult(temp,temp,temp);
658: VecDot(temp,appctx->SEMop.mass,&gnorm);
659: gnorm = PetscSqrtReal(gnorm);
660: VecDestroy(&temp);
661: TaoGetIterationNumber(tao,&its);
662: TaoGetObjective(tao,&fct);
663: if (!its) {
664: PetscPrintf(PETSC_COMM_WORLD,"%% Iteration Error Objective Gradient-norm\n");
665: PetscPrintf(PETSC_COMM_WORLD,"history = [\n");
666: }
667: PetscPrintf(PETSC_COMM_WORLD,"%3D %g %g %g\n",its,(double)nrm,(double)fct,(double)gnorm);
668: return(0);
669: }
671: PetscErrorCode MonitorDestroy(void **ctx)
672: {
676: PetscPrintf(PETSC_COMM_WORLD,"];\n");
677: return(0);
678: }
680: /*TEST
682: build:
683: requires: !complex
685: test:
686: requires: !single
687: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
689: test:
690: suffix: cn
691: requires: !single
692: args: -ts_type cn -ts_dt .003 -pc_type lu -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
694: test:
695: suffix: 2
696: requires: !single
697: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -a .1 -tao_bqnls_mat_lmvm_scale_type none
699: TEST*/