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 */
119: PetscInt i, xs, xm, ind, j, lenglob;
120: PetscReal x, *wrk_ptr1, *wrk_ptr2;
121: MatNullSpace nsp;
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: Initialize program and set problem parameters
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127: PetscInitialize(&argc,&argv,(char*)0,help);
129: /*initialize parameters */
130: appctx.param.N = 10; /* order of the spectral element */
131: appctx.param.E = 8; /* number of elements */
132: appctx.param.L = 1.0; /* length of the domain */
133: appctx.param.mu = 0.00001; /* diffusion coefficient */
134: appctx.param.a = 0.0; /* advection speed */
135: appctx.initial_dt = 1e-4;
136: appctx.param.steps = PETSC_MAX_INT;
137: appctx.param.Tend = 0.01;
138: appctx.ncoeff = 2;
140: PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL);
141: PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL);
142: PetscOptionsGetInt(NULL,NULL,"-ncoeff",&appctx.ncoeff,NULL);
143: PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL);
144: PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL);
145: PetscOptionsGetReal(NULL,NULL,"-a",&appctx.param.a,NULL);
146: appctx.param.Le = appctx.param.L/appctx.param.E;
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Create GLL data structures
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: PetscMalloc2(appctx.param.N,&appctx.SEMop.gll.nodes,appctx.param.N,&appctx.SEMop.gll.weights);
152: PetscDTGaussLobattoLegendreQuadrature(appctx.param.N,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
153: appctx.SEMop.gll.n = appctx.param.N;
154: lenglob = appctx.param.E*(appctx.param.N-1);
156: /*
157: Create distributed array (DMDA) to manage parallel grid and vectors
158: and to set up the ghost point communication pattern. There are E*(Nl-1)+1
159: total grid values spread equally among all the processors, except first and last
160: */
162: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);
163: DMSetFromOptions(appctx.da);
164: DMSetUp(appctx.da);
166: /*
167: Extract global and local vectors from DMDA; we use these to store the
168: approximate solution. Then duplicate these for remaining vectors that
169: have the same types.
170: */
172: DMCreateGlobalVector(appctx.da,&u);
173: VecDuplicate(u,&appctx.dat.ic);
174: VecDuplicate(u,&appctx.dat.true_solution);
175: VecDuplicate(u,&appctx.dat.reference);
176: VecDuplicate(u,&appctx.SEMop.grid);
177: VecDuplicate(u,&appctx.SEMop.mass);
178: VecDuplicate(u,&appctx.dat.curr_sol);
179: VecDuplicate(u,&appctx.dat.joe);
181: DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
182: DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
183: DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
185: /* Compute function over the locally owned part of the grid */
187: xs=xs/(appctx.param.N-1);
188: xm=xm/(appctx.param.N-1);
190: /*
191: Build total grid and mass over entire mesh (multi-elemental)
192: */
194: for (i=xs; i<xs+xm; i++) {
195: for (j=0; j<appctx.param.N-1; j++) {
196: x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i;
197: ind=i*(appctx.param.N-1)+j;
198: wrk_ptr1[ind]=x;
199: wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
200: if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
201: }
202: }
203: DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
204: DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
206: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207: Create matrix data structure; set matrix evaluation routine.
208: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209: DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);
210: DMCreateMatrix(appctx.da,&appctx.SEMop.stiff);
211: DMCreateMatrix(appctx.da,&appctx.SEMop.advec);
213: /*
214: For linear problems with a time-dependent f(u,t) in the equation
215: u_t = f(u,t), the user provides the discretized right-hand-side
216: as a time-dependent matrix.
217: */
218: RHSLaplacian(appctx.ts,0.0,u,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx);
219: RHSAdvection(appctx.ts,0.0,u,appctx.SEMop.advec,appctx.SEMop.advec,&appctx);
220: MatAXPY(appctx.SEMop.stiff,-1.0,appctx.SEMop.advec,DIFFERENT_NONZERO_PATTERN);
221: MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff);
223: /* attach the null space to the matrix, this probably is not needed but does no harm */
224: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
225: MatSetNullSpace(appctx.SEMop.stiff,nsp);
226: MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL);
227: MatNullSpaceDestroy(&nsp);
229: /* Create the TS solver that solves the ODE and its adjoint; set its options */
230: TSCreate(PETSC_COMM_WORLD,&appctx.ts);
231: TSSetSolutionFunction(appctx.ts,(PetscErrorCode (*)(TS,PetscReal,Vec, void *))ComputeReference,&appctx);
232: TSSetProblemType(appctx.ts,TS_LINEAR);
233: TSSetType(appctx.ts,TSRK);
234: TSSetDM(appctx.ts,appctx.da);
235: TSSetTime(appctx.ts,0.0);
236: TSSetTimeStep(appctx.ts,appctx.initial_dt);
237: TSSetMaxSteps(appctx.ts,appctx.param.steps);
238: TSSetMaxTime(appctx.ts,appctx.param.Tend);
239: TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
240: TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL);
241: TSSetFromOptions(appctx.ts);
242: /* Need to save initial timestep user may have set with -ts_dt so it can be reset for each new TSSolve() */
243: TSGetTimeStep(appctx.ts,&appctx.initial_dt);
244: TSSetRHSFunction(appctx.ts,NULL,TSComputeRHSFunctionLinear,&appctx);
245: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,TSComputeRHSJacobianConstant,&appctx);
246: /* TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
247: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx); */
249: /* Set random initial conditions as initial guess, compute analytic reference solution and analytic (true) initial conditions */
250: ComputeSolutionCoefficients(&appctx);
251: InitialConditions(appctx.dat.ic,&appctx);
252: ComputeReference(appctx.ts,appctx.param.Tend,appctx.dat.reference,&appctx);
253: ComputeReference(appctx.ts,0.0,appctx.dat.true_solution,&appctx);
255: /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */
256: TSSetSaveTrajectory(appctx.ts);
257: TSSetFromOptions(appctx.ts);
259: /* Create TAO solver and set desired solution method */
260: TaoCreate(PETSC_COMM_WORLD,&tao);
261: TaoSetMonitor(tao,MonitorError,&appctx,MonitorDestroy);
262: TaoSetType(tao,TAOBQNLS);
263: TaoSetSolution(tao,appctx.dat.ic);
264: /* Set routine for function and gradient evaluation */
265: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&appctx);
266: /* Check for any TAO command line options */
267: TaoSetTolerances(tao,1e-8,PETSC_DEFAULT,PETSC_DEFAULT);
268: TaoSetFromOptions(tao);
269: TaoSolve(tao);
271: TaoDestroy(&tao);
272: PetscFree(appctx.solutioncoefficients);
273: MatDestroy(&appctx.SEMop.advec);
274: MatDestroy(&appctx.SEMop.stiff);
275: MatDestroy(&appctx.SEMop.keptstiff);
276: VecDestroy(&u);
277: VecDestroy(&appctx.dat.ic);
278: VecDestroy(&appctx.dat.joe);
279: VecDestroy(&appctx.dat.true_solution);
280: VecDestroy(&appctx.dat.reference);
281: VecDestroy(&appctx.SEMop.grid);
282: VecDestroy(&appctx.SEMop.mass);
283: VecDestroy(&appctx.dat.curr_sol);
284: PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
285: DMDestroy(&appctx.da);
286: TSDestroy(&appctx.ts);
288: /*
289: Always call PetscFinalize() before exiting a program. This routine
290: - finalizes the PETSc libraries as well as MPI
291: - provides summary and diagnostic information if certain runtime
292: options are chosen (e.g., -log_summary).
293: */
294: PetscFinalize();
295: return 0;
296: }
298: /*
299: Computes the coefficients for the analytic solution to the PDE
300: */
301: PetscErrorCode ComputeSolutionCoefficients(AppCtx *appctx)
302: {
303: PetscRandom rand;
304: PetscInt i;
306: PetscMalloc1(appctx->ncoeff,&appctx->solutioncoefficients);
307: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
308: PetscRandomSetInterval(rand,.9,1.0);
309: for (i=0; i<appctx->ncoeff; i++) {
310: PetscRandomGetValue(rand,&appctx->solutioncoefficients[i]);
311: }
312: PetscRandomDestroy(&rand);
313: return 0;
314: }
316: /* --------------------------------------------------------------------- */
317: /*
318: InitialConditions - Computes the (random) initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()
320: Input Parameter:
321: u - uninitialized solution vector (global)
322: appctx - user-defined application context
324: Output Parameter:
325: u - vector with solution at initial time (global)
326: */
327: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
328: {
329: PetscScalar *s;
330: const PetscScalar *xg;
331: PetscInt i,j,lenglob;
332: PetscReal sum,val;
333: PetscRandom rand;
335: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
336: PetscRandomSetInterval(rand,.9,1.0);
337: DMDAVecGetArray(appctx->da,u,&s);
338: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
339: lenglob = appctx->param.E*(appctx->param.N-1);
340: for (i=0; i<lenglob; i++) {
341: s[i]= 0;
342: for (j=0; j<appctx->ncoeff; j++) {
343: PetscRandomGetValue(rand,&val);
344: s[i] += val*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
345: }
346: }
347: PetscRandomDestroy(&rand);
348: DMDAVecRestoreArray(appctx->da,u,&s);
349: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
350: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
351: VecSum(u,&sum);
352: VecShift(u,-sum/lenglob);
353: return 0;
354: }
356: /*
357: TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function.
359: InitialConditions() computes the initial conditions for the beginning of the Tao iterations
361: Input Parameter:
362: u - uninitialized solution vector (global)
363: appctx - user-defined application context
365: Output Parameter:
366: u - vector with solution at initial time (global)
367: */
368: PetscErrorCode TrueSolution(Vec u,AppCtx *appctx)
369: {
370: PetscScalar *s;
371: const PetscScalar *xg;
372: PetscInt i,j,lenglob;
373: 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: PetscInt i, j,lenglob;
407: DMDAVecGetArray(appctx->da,obj,&s);
408: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
409: lenglob = appctx->param.E*(appctx->param.N-1);
410: for (i=0; i<lenglob; i++) {
411: s[i]= 0;
412: for (j=0; j<appctx->ncoeff; j++) {
413: tc = -appctx->param.mu*(j+1)*(j+1)*4.0*PETSC_PI*PETSC_PI*t;
414: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*(xg[i] + appctx->param.a*t))*PetscExpReal(tc);
415: }
416: }
417: DMDAVecRestoreArray(appctx->da,obj,&s);
418: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
419: return 0;
420: }
422: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
423: {
424: AppCtx *appctx = (AppCtx*)ctx;
426: MatMult(appctx->SEMop.keptstiff,globalin,globalout);
427: return 0;
428: }
430: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx)
431: {
432: AppCtx *appctx = (AppCtx*)ctx;
434: MatCopy(appctx->SEMop.keptstiff,A,DIFFERENT_NONZERO_PATTERN);
435: return 0;
436: }
438: /* --------------------------------------------------------------------- */
440: /*
441: RHSLaplacian - matrix for diffusion
443: Input Parameters:
444: ts - the TS context
445: t - current time (ignored)
446: X - current solution (ignored)
447: dummy - optional user-defined context, as set by TSetRHSJacobian()
449: Output Parameters:
450: AA - Jacobian matrix
451: BB - optionally different matrix from which the preconditioner is built
452: str - flag indicating matrix structure
454: Scales by the inverse of the mass matrix (perhaps that should be pulled out)
456: */
457: PetscErrorCode RHSLaplacian(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
458: {
459: PetscReal **temp;
460: PetscReal vv;
461: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
462: PetscInt i,xs,xn,l,j;
463: PetscInt *rowsDM;
465: /*
466: Creates the element stiffness matrix for the given gll
467: */
468: PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
470: /* scale by the size of the element */
471: for (i=0; i<appctx->param.N; i++) {
472: vv=-appctx->param.mu*2.0/appctx->param.Le;
473: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
474: }
476: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
477: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
480: xs = xs/(appctx->param.N-1);
481: xn = xn/(appctx->param.N-1);
483: PetscMalloc1(appctx->param.N,&rowsDM);
484: /*
485: loop over local elements
486: */
487: for (j=xs; j<xs+xn; j++) {
488: for (l=0; l<appctx->param.N; l++) {
489: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
490: }
491: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
492: }
493: PetscFree(rowsDM);
494: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
495: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
496: VecReciprocal(appctx->SEMop.mass);
497: MatDiagonalScale(A,appctx->SEMop.mass,0);
498: VecReciprocal(appctx->SEMop.mass);
500: PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
501: return 0;
502: }
504: /*
505: Almost identical to Laplacian
507: Note that the element matrix is NOT scaled by the size of element like the Laplacian term.
508: */
509: PetscErrorCode RHSAdvection(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
510: {
511: PetscReal **temp;
512: PetscReal vv;
513: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
514: PetscInt i,xs,xn,l,j;
515: PetscInt *rowsDM;
517: /*
518: Creates the element stiffness matrix for the given gll
519: */
520: PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
522: /* scale by the size of the element */
523: for (i=0; i<appctx->param.N; i++) {
524: vv = -appctx->param.a;
525: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
526: }
528: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
529: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
532: xs = xs/(appctx->param.N-1);
533: xn = xn/(appctx->param.N-1);
535: PetscMalloc1(appctx->param.N,&rowsDM);
536: /*
537: loop over local elements
538: */
539: for (j=xs; j<xs+xn; j++) {
540: for (l=0; l<appctx->param.N; l++) {
541: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
542: }
543: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
544: }
545: PetscFree(rowsDM);
546: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
547: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
548: VecReciprocal(appctx->SEMop.mass);
549: MatDiagonalScale(A,appctx->SEMop.mass,0);
550: VecReciprocal(appctx->SEMop.mass);
552: PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
553: return 0;
554: }
556: /* ------------------------------------------------------------------ */
557: /*
558: FormFunctionGradient - Evaluates the function and corresponding gradient.
560: Input Parameters:
561: tao - the Tao context
562: ic - the input vector
563: ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradient()
565: Output Parameters:
566: f - the newly evaluated function
567: G - the newly evaluated gradient
569: Notes:
571: The forward equation is
572: M u_t = F(U)
573: which is converted to
574: u_t = M^{-1} F(u)
575: in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is
576: M^{-1} J
577: where J is the Jacobian of F. Now the adjoint equation is
578: M v_t = J^T v
579: but TSAdjoint does not solve this since it can only solve the transposed system for the
580: Jacobian the user provided. Hence TSAdjoint solves
581: w_t = J^T M^{-1} w (where w = M v)
582: since there is no way to indicate the mass matrix as a separate entity to TS. Thus one
583: must be careful in initializing the "adjoint equation" and using the result. This is
584: why
585: G = -2 M(u(T) - u_d)
586: below (instead of -2(u(T) - u_d)
588: */
589: PetscErrorCode FormFunctionGradient(Tao tao,Vec ic,PetscReal *f,Vec G,void *ctx)
590: {
591: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
592: Vec temp;
594: TSSetTime(appctx->ts,0.0);
595: TSSetStepNumber(appctx->ts,0);
596: TSSetTimeStep(appctx->ts,appctx->initial_dt);
597: VecCopy(ic,appctx->dat.curr_sol);
599: TSSolve(appctx->ts,appctx->dat.curr_sol);
600: VecCopy(appctx->dat.curr_sol,appctx->dat.joe);
602: /* Compute the difference between the current ODE solution and target ODE solution */
603: VecWAXPY(G,-1.0,appctx->dat.curr_sol,appctx->dat.reference);
605: /* Compute the objective/cost function */
606: VecDuplicate(G,&temp);
607: VecPointwiseMult(temp,G,G);
608: VecDot(temp,appctx->SEMop.mass,f);
609: VecDestroy(&temp);
611: /* Compute initial conditions for the adjoint integration. See Notes above */
612: VecScale(G, -2.0);
613: VecPointwiseMult(G,G,appctx->SEMop.mass);
614: TSSetCostGradients(appctx->ts,1,&G,NULL);
616: TSAdjointSolve(appctx->ts);
617: /* VecPointwiseDivide(G,G,appctx->SEMop.mass);*/
618: return 0;
619: }
621: PetscErrorCode MonitorError(Tao tao,void *ctx)
622: {
623: AppCtx *appctx = (AppCtx*)ctx;
624: Vec temp,grad;
625: PetscReal nrm;
626: PetscInt its;
627: PetscReal fct,gnorm;
629: VecDuplicate(appctx->dat.ic,&temp);
630: VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution);
631: VecPointwiseMult(temp,temp,temp);
632: VecDot(temp,appctx->SEMop.mass,&nrm);
633: nrm = PetscSqrtReal(nrm);
634: TaoGetGradient(tao,&grad,NULL,NULL);
635: VecPointwiseMult(temp,temp,temp);
636: VecDot(temp,appctx->SEMop.mass,&gnorm);
637: gnorm = PetscSqrtReal(gnorm);
638: VecDestroy(&temp);
639: TaoGetIterationNumber(tao,&its);
640: TaoGetSolutionStatus(tao,NULL,&fct,NULL,NULL,NULL,NULL);
641: if (!its) {
642: PetscPrintf(PETSC_COMM_WORLD,"%% Iteration Error Objective Gradient-norm\n");
643: PetscPrintf(PETSC_COMM_WORLD,"history = [\n");
644: }
645: PetscPrintf(PETSC_COMM_WORLD,"%3D %g %g %g\n",its,(double)nrm,(double)fct,(double)gnorm);
646: return 0;
647: }
649: PetscErrorCode MonitorDestroy(void **ctx)
650: {
651: PetscPrintf(PETSC_COMM_WORLD,"];\n");
652: return 0;
653: }
655: /*TEST
657: build:
658: requires: !complex
660: test:
661: requires: !single
662: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
664: test:
665: suffix: cn
666: requires: !single
667: args: -ts_type cn -ts_dt .003 -pc_type lu -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
669: test:
670: suffix: 2
671: requires: !single
672: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -a .1 -tao_bqnls_mat_lmvm_scale_type none
674: TEST*/