Actual source code: spectraladjointassimilation.c
petsc-3.12.5 2020-03-29
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;
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Create GLL data structures
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: PetscMalloc2(appctx.param.N,&appctx.SEMop.gll.nodes,appctx.param.N,&appctx.SEMop.gll.weights);
155: PetscDTGaussLobattoLegendreQuadrature(appctx.param.N,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
156: appctx.SEMop.gll.n = appctx.param.N;
157: lenglob = appctx.param.E*(appctx.param.N-1);
159: /*
160: Create distributed array (DMDA) to manage parallel grid and vectors
161: and to set up the ghost point communication pattern. There are E*(Nl-1)+1
162: total grid values spread equally among all the processors, except first and last
163: */
165: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);
166: DMSetFromOptions(appctx.da);
167: DMSetUp(appctx.da);
169: /*
170: Extract global and local vectors from DMDA; we use these to store the
171: approximate solution. Then duplicate these for remaining vectors that
172: have the same types.
173: */
175: DMCreateGlobalVector(appctx.da,&u);
176: VecDuplicate(u,&appctx.dat.ic);
177: VecDuplicate(u,&appctx.dat.true_solution);
178: VecDuplicate(u,&appctx.dat.reference);
179: VecDuplicate(u,&appctx.SEMop.grid);
180: VecDuplicate(u,&appctx.SEMop.mass);
181: VecDuplicate(u,&appctx.dat.curr_sol);
182: VecDuplicate(u,&appctx.dat.joe);
184: DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
185: DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
186: DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
188: /* Compute function over the locally owned part of the grid */
190: xs=xs/(appctx.param.N-1);
191: xm=xm/(appctx.param.N-1);
193: /*
194: Build total grid and mass over entire mesh (multi-elemental)
195: */
197: for (i=xs; i<xs+xm; i++) {
198: for (j=0; j<appctx.param.N-1; j++) {
199: x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i;
200: ind=i*(appctx.param.N-1)+j;
201: wrk_ptr1[ind]=x;
202: wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
203: if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
204: }
205: }
206: DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
207: DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
210: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211: Create matrix data structure; set matrix evaluation routine.
212: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
213: DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);
214: DMCreateMatrix(appctx.da,&appctx.SEMop.stiff);
215: DMCreateMatrix(appctx.da,&appctx.SEMop.advec);
217: /*
218: For linear problems with a time-dependent f(u,t) in the equation
219: u_t = f(u,t), the user provides the discretized right-hand-side
220: as a time-dependent matrix.
221: */
222: RHSLaplacian(appctx.ts,0.0,u,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx);
223: RHSAdvection(appctx.ts,0.0,u,appctx.SEMop.advec,appctx.SEMop.advec,&appctx);
224: MatAXPY(appctx.SEMop.stiff,-1.0,appctx.SEMop.advec,DIFFERENT_NONZERO_PATTERN);
225: MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff);
227: /* attach the null space to the matrix, this probably is not needed but does no harm */
228: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
229: MatSetNullSpace(appctx.SEMop.stiff,nsp);
230: MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL);
231: MatNullSpaceDestroy(&nsp);
233: /* Create the TS solver that solves the ODE and its adjoint; set its options */
234: TSCreate(PETSC_COMM_WORLD,&appctx.ts);
235: TSSetSolutionFunction(appctx.ts,(PetscErrorCode (*)(TS,PetscReal,Vec, void *))ComputeReference,&appctx);
236: TSSetProblemType(appctx.ts,TS_LINEAR);
237: TSSetType(appctx.ts,TSRK);
238: TSSetDM(appctx.ts,appctx.da);
239: TSSetTime(appctx.ts,0.0);
240: TSSetTimeStep(appctx.ts,appctx.initial_dt);
241: TSSetMaxSteps(appctx.ts,appctx.param.steps);
242: TSSetMaxTime(appctx.ts,appctx.param.Tend);
243: TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
244: TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL);
245: TSSetFromOptions(appctx.ts);
246: /* Need to save initial timestep user may have set with -ts_dt so it can be reset for each new TSSolve() */
247: TSGetTimeStep(appctx.ts,&appctx.initial_dt);
248: TSSetRHSFunction(appctx.ts,NULL,TSComputeRHSFunctionLinear,&appctx);
249: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,TSComputeRHSJacobianConstant,&appctx);
250: /* TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
251: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx); */
252: TSSetSaveTrajectory(appctx.ts);
254: /* Set random initial conditions as initial guess, compute analytic reference solution and analytic (true) initial conditions */
255: ComputeSolutionCoefficients(&appctx);
256: InitialConditions(appctx.dat.ic,&appctx);
257: ComputeReference(appctx.ts,appctx.param.Tend,appctx.dat.reference,&appctx);
258: ComputeReference(appctx.ts,0.0,appctx.dat.true_solution,&appctx);
260: /* Create TAO solver and set desired solution method */
261: TaoCreate(PETSC_COMM_WORLD,&tao);
262: TaoSetMonitor(tao,MonitorError,&appctx,MonitorDestroy);
263: TaoSetType(tao,TAOBQNLS);
264: TaoSetInitialVector(tao,appctx.dat.ic);
265: /* Set routine for function and gradient evaluation */
266: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&appctx);
267: /* Check for any TAO command line options */
268: TaoSetTolerances(tao,1e-8,PETSC_DEFAULT,PETSC_DEFAULT);
269: TaoSetFromOptions(tao);
270: TaoSolve(tao);
272: TaoDestroy(&tao);
273: PetscFree(appctx.solutioncoefficients);
274: MatDestroy(&appctx.SEMop.advec);
275: MatDestroy(&appctx.SEMop.stiff);
276: MatDestroy(&appctx.SEMop.keptstiff);
277: VecDestroy(&u);
278: VecDestroy(&appctx.dat.ic);
279: VecDestroy(&appctx.dat.joe);
280: VecDestroy(&appctx.dat.true_solution);
281: VecDestroy(&appctx.dat.reference);
282: VecDestroy(&appctx.SEMop.grid);
283: VecDestroy(&appctx.SEMop.mass);
284: VecDestroy(&appctx.dat.curr_sol);
285: PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
286: DMDestroy(&appctx.da);
287: TSDestroy(&appctx.ts);
289: /*
290: Always call PetscFinalize() before exiting a program. This routine
291: - finalizes the PETSc libraries as well as MPI
292: - provides summary and diagnostic information if certain runtime
293: options are chosen (e.g., -log_summary).
294: */
295: PetscFinalize();
296: return ierr;
297: }
299: /*
300: Computes the coefficients for the analytic solution to the PDE
301: */
302: PetscErrorCode ComputeSolutionCoefficients(AppCtx *appctx)
303: {
304: PetscErrorCode ierr;
305: PetscRandom rand;
306: PetscInt i;
309: PetscMalloc1(appctx->ncoeff,&appctx->solutioncoefficients);
310: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
311: PetscRandomSetInterval(rand,.9,1.0);
312: for (i=0; i<appctx->ncoeff; i++) {
313: PetscRandomGetValue(rand,&appctx->solutioncoefficients[i]);
314: }
315: PetscRandomDestroy(&rand);
316: return(0);
317: }
319: /* --------------------------------------------------------------------- */
320: /*
321: InitialConditions - Computes the (random) initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()
323: Input Parameter:
324: u - uninitialized solution vector (global)
325: appctx - user-defined application context
327: Output Parameter:
328: u - vector with solution at initial time (global)
329: */
330: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
331: {
332: PetscScalar *s;
333: const PetscScalar *xg;
334: PetscErrorCode ierr;
335: PetscInt i,j,lenglob;
336: PetscReal sum,val;
337: PetscRandom rand;
340: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
341: PetscRandomSetInterval(rand,.9,1.0);
342: DMDAVecGetArray(appctx->da,u,&s);
343: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
344: lenglob = appctx->param.E*(appctx->param.N-1);
345: for (i=0; i<lenglob; i++) {
346: s[i]= 0;
347: for (j=0; j<appctx->ncoeff; j++) {
348: PetscRandomGetValue(rand,&val);
349: s[i] += val*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
350: }
351: }
352: PetscRandomDestroy(&rand);
353: DMDAVecRestoreArray(appctx->da,u,&s);
354: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
355: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
356: VecSum(u,&sum);
357: VecShift(u,-sum/lenglob);
358: return(0);
359: }
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 begining 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 entitity 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)
607: */
608: PetscErrorCode FormFunctionGradient(Tao tao,Vec ic,PetscReal *f,Vec G,void *ctx)
609: {
610: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
611: PetscErrorCode ierr;
612: Vec temp;
615: TSSetTime(appctx->ts,0.0);
616: TSSetStepNumber(appctx->ts,0);
617: TSResetTrajectory(appctx->ts);
618: TSSetTimeStep(appctx->ts,appctx->initial_dt);
619: VecCopy(ic,appctx->dat.curr_sol);
621: TSSolve(appctx->ts,appctx->dat.curr_sol);
622: VecCopy(appctx->dat.curr_sol,appctx->dat.joe);
624: /* Compute the difference between the current ODE solution and target ODE solution */
625: VecWAXPY(G,-1.0,appctx->dat.curr_sol,appctx->dat.reference);
627: /* Compute the objective/cost function */
628: VecDuplicate(G,&temp);
629: VecPointwiseMult(temp,G,G);
630: VecDot(temp,appctx->SEMop.mass,f);
631: VecDestroy(&temp);
633: /* Compute initial conditions for the adjoint integration. See Notes above */
634: VecScale(G, -2.0);
635: VecPointwiseMult(G,G,appctx->SEMop.mass);
636: TSSetCostGradients(appctx->ts,1,&G,NULL);
638: TSAdjointSolve(appctx->ts);
639: /* VecPointwiseDivide(G,G,appctx->SEMop.mass);*/
640: return(0);
641: }
643: PetscErrorCode MonitorError(Tao tao,void *ctx)
644: {
645: AppCtx *appctx = (AppCtx*)ctx;
646: Vec temp,grad;
647: PetscReal nrm;
649: PetscInt its;
650: PetscReal fct,gnorm;
653: VecDuplicate(appctx->dat.ic,&temp);
654: VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution);
655: VecPointwiseMult(temp,temp,temp);
656: VecDot(temp,appctx->SEMop.mass,&nrm);
657: nrm = PetscSqrtReal(nrm);
658: TaoGetGradientVector(tao,&grad);
659: VecPointwiseMult(temp,temp,temp);
660: VecDot(temp,appctx->SEMop.mass,&gnorm);
661: gnorm = PetscSqrtReal(gnorm);
662: VecDestroy(&temp);
663: TaoGetIterationNumber(tao,&its);
664: TaoGetObjective(tao,&fct);
665: if (!its) {
666: PetscPrintf(PETSC_COMM_WORLD,"%% Iteration Error Objective Gradient-norm\n");
667: PetscPrintf(PETSC_COMM_WORLD,"history = [\n");
668: }
669: PetscPrintf(PETSC_COMM_WORLD,"%3D %g %g %g\n",its,(double)nrm,(double)fct,(double)gnorm);
670: return(0);
671: }
673: PetscErrorCode MonitorDestroy(void **ctx)
674: {
678: PetscPrintf(PETSC_COMM_WORLD,"];\n");
679: return(0);
680: }
683: /*TEST
685: build:
686: requires: !complex
688: test:
689: requires: !single
690: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
692: test:
693: suffix: cn
694: requires: !single
695: args: -ts_type cn -ts_dt .003 -pc_type lu -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
697: test:
698: suffix: 2
699: requires: !single
700: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -a .1 -tao_bqnls_mat_lmvm_scale_type none
703: TEST*/