Actual source code: spectraladjointassimilation.c
petsc-3.14.6 2021-03-30
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); */
253: /* Set random initial conditions as initial guess, compute analytic reference solution and analytic (true) initial conditions */
254: ComputeSolutionCoefficients(&appctx);
255: InitialConditions(appctx.dat.ic,&appctx);
256: ComputeReference(appctx.ts,appctx.param.Tend,appctx.dat.reference,&appctx);
257: ComputeReference(appctx.ts,0.0,appctx.dat.true_solution,&appctx);
259: /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */
260: TSSetSaveTrajectory(appctx.ts);
261: TSSetFromOptions(appctx.ts);
263: /* Create TAO solver and set desired solution method */
264: TaoCreate(PETSC_COMM_WORLD,&tao);
265: TaoSetMonitor(tao,MonitorError,&appctx,MonitorDestroy);
266: TaoSetType(tao,TAOBQNLS);
267: TaoSetInitialVector(tao,appctx.dat.ic);
268: /* Set routine for function and gradient evaluation */
269: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&appctx);
270: /* Check for any TAO command line options */
271: TaoSetTolerances(tao,1e-8,PETSC_DEFAULT,PETSC_DEFAULT);
272: TaoSetFromOptions(tao);
273: TaoSolve(tao);
275: TaoDestroy(&tao);
276: PetscFree(appctx.solutioncoefficients);
277: MatDestroy(&appctx.SEMop.advec);
278: MatDestroy(&appctx.SEMop.stiff);
279: MatDestroy(&appctx.SEMop.keptstiff);
280: VecDestroy(&u);
281: VecDestroy(&appctx.dat.ic);
282: VecDestroy(&appctx.dat.joe);
283: VecDestroy(&appctx.dat.true_solution);
284: VecDestroy(&appctx.dat.reference);
285: VecDestroy(&appctx.SEMop.grid);
286: VecDestroy(&appctx.SEMop.mass);
287: VecDestroy(&appctx.dat.curr_sol);
288: PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
289: DMDestroy(&appctx.da);
290: TSDestroy(&appctx.ts);
292: /*
293: Always call PetscFinalize() before exiting a program. This routine
294: - finalizes the PETSc libraries as well as MPI
295: - provides summary and diagnostic information if certain runtime
296: options are chosen (e.g., -log_summary).
297: */
298: PetscFinalize();
299: return ierr;
300: }
302: /*
303: Computes the coefficients for the analytic solution to the PDE
304: */
305: PetscErrorCode ComputeSolutionCoefficients(AppCtx *appctx)
306: {
307: PetscErrorCode ierr;
308: PetscRandom rand;
309: PetscInt i;
312: PetscMalloc1(appctx->ncoeff,&appctx->solutioncoefficients);
313: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
314: PetscRandomSetInterval(rand,.9,1.0);
315: for (i=0; i<appctx->ncoeff; i++) {
316: PetscRandomGetValue(rand,&appctx->solutioncoefficients[i]);
317: }
318: PetscRandomDestroy(&rand);
319: return(0);
320: }
322: /* --------------------------------------------------------------------- */
323: /*
324: InitialConditions - Computes the (random) initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()
326: Input Parameter:
327: u - uninitialized solution vector (global)
328: appctx - user-defined application context
330: Output Parameter:
331: u - vector with solution at initial time (global)
332: */
333: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
334: {
335: PetscScalar *s;
336: const PetscScalar *xg;
337: PetscErrorCode ierr;
338: PetscInt i,j,lenglob;
339: PetscReal sum,val;
340: PetscRandom rand;
343: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
344: PetscRandomSetInterval(rand,.9,1.0);
345: DMDAVecGetArray(appctx->da,u,&s);
346: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
347: lenglob = appctx->param.E*(appctx->param.N-1);
348: for (i=0; i<lenglob; i++) {
349: s[i]= 0;
350: for (j=0; j<appctx->ncoeff; j++) {
351: PetscRandomGetValue(rand,&val);
352: s[i] += val*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
353: }
354: }
355: PetscRandomDestroy(&rand);
356: DMDAVecRestoreArray(appctx->da,u,&s);
357: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
358: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
359: VecSum(u,&sum);
360: VecShift(u,-sum/lenglob);
361: return(0);
362: }
365: /*
366: TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function.
368: InitialConditions() computes the initial conditions for the begining of the Tao iterations
370: Input Parameter:
371: u - uninitialized solution vector (global)
372: appctx - user-defined application context
374: Output Parameter:
375: u - vector with solution at initial time (global)
376: */
377: PetscErrorCode TrueSolution(Vec u,AppCtx *appctx)
378: {
379: PetscScalar *s;
380: const PetscScalar *xg;
381: PetscErrorCode ierr;
382: PetscInt i,j,lenglob;
383: PetscReal sum;
386: DMDAVecGetArray(appctx->da,u,&s);
387: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
388: lenglob = appctx->param.E*(appctx->param.N-1);
389: for (i=0; i<lenglob; i++) {
390: s[i]= 0;
391: for (j=0; j<appctx->ncoeff; j++) {
392: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
393: }
394: }
395: DMDAVecRestoreArray(appctx->da,u,&s);
396: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
397: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
398: VecSum(u,&sum);
399: VecShift(u,-sum/lenglob);
400: return(0);
401: }
402: /* --------------------------------------------------------------------- */
403: /*
404: Sets the desired profile for the final end time
406: Input Parameters:
407: t - final time
408: obj - vector storing the desired profile
409: appctx - user-defined application context
411: */
412: PetscErrorCode ComputeReference(TS ts,PetscReal t,Vec obj,AppCtx *appctx)
413: {
414: PetscScalar *s,tc;
415: const PetscScalar *xg;
416: PetscErrorCode ierr;
417: PetscInt i, j,lenglob;
420: DMDAVecGetArray(appctx->da,obj,&s);
421: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
422: lenglob = appctx->param.E*(appctx->param.N-1);
423: for (i=0; i<lenglob; i++) {
424: s[i]= 0;
425: for (j=0; j<appctx->ncoeff; j++) {
426: tc = -appctx->param.mu*(j+1)*(j+1)*4.0*PETSC_PI*PETSC_PI*t;
427: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*(xg[i] + appctx->param.a*t))*PetscExpReal(tc);
428: }
429: }
430: DMDAVecRestoreArray(appctx->da,obj,&s);
431: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
432: return(0);
433: }
435: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
436: {
438: AppCtx *appctx = (AppCtx*)ctx;
441: MatMult(appctx->SEMop.keptstiff,globalin,globalout);
442: return(0);
443: }
445: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx)
446: {
448: AppCtx *appctx = (AppCtx*)ctx;
451: MatCopy(appctx->SEMop.keptstiff,A,DIFFERENT_NONZERO_PATTERN);
452: return(0);
453: }
455: /* --------------------------------------------------------------------- */
457: /*
458: RHSLaplacian - matrix for diffusion
460: Input Parameters:
461: ts - the TS context
462: t - current time (ignored)
463: X - current solution (ignored)
464: dummy - optional user-defined context, as set by TSetRHSJacobian()
466: Output Parameters:
467: AA - Jacobian matrix
468: BB - optionally different matrix from which the preconditioner is built
469: str - flag indicating matrix structure
471: Scales by the inverse of the mass matrix (perhaps that should be pulled out)
473: */
474: PetscErrorCode RHSLaplacian(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
475: {
476: PetscReal **temp;
477: PetscReal vv;
478: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
480: PetscInt i,xs,xn,l,j;
481: PetscInt *rowsDM;
484: /*
485: Creates the element stiffness matrix for the given gll
486: */
487: PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
489: /* scale by the size of the element */
490: for (i=0; i<appctx->param.N; i++) {
491: vv=-appctx->param.mu*2.0/appctx->param.Le;
492: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
493: }
495: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
496: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
498: if (appctx->param.N-1 < 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Polynomial order must be at least 2");
499: xs = xs/(appctx->param.N-1);
500: xn = xn/(appctx->param.N-1);
502: PetscMalloc1(appctx->param.N,&rowsDM);
503: /*
504: loop over local elements
505: */
506: for (j=xs; j<xs+xn; j++) {
507: for (l=0; l<appctx->param.N; l++) {
508: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
509: }
510: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
511: }
512: PetscFree(rowsDM);
513: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
514: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
515: VecReciprocal(appctx->SEMop.mass);
516: MatDiagonalScale(A,appctx->SEMop.mass,0);
517: VecReciprocal(appctx->SEMop.mass);
519: PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
520: return(0);
521: }
523: /*
524: Almost identical to Laplacian
526: Note that the element matrix is NOT scaled by the size of element like the Laplacian term.
527: */
528: PetscErrorCode RHSAdvection(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
529: {
530: PetscReal **temp;
531: PetscReal vv;
532: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
534: PetscInt i,xs,xn,l,j;
535: PetscInt *rowsDM;
538: /*
539: Creates the element stiffness matrix for the given gll
540: */
541: PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
543: /* scale by the size of the element */
544: for (i=0; i<appctx->param.N; i++) {
545: vv = -appctx->param.a;
546: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
547: }
549: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
550: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
552: if (appctx->param.N-1 < 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Polynomial order must be at least 2");
553: xs = xs/(appctx->param.N-1);
554: xn = xn/(appctx->param.N-1);
556: PetscMalloc1(appctx->param.N,&rowsDM);
557: /*
558: loop over local elements
559: */
560: for (j=xs; j<xs+xn; j++) {
561: for (l=0; l<appctx->param.N; l++) {
562: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
563: }
564: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
565: }
566: PetscFree(rowsDM);
567: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
568: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
569: VecReciprocal(appctx->SEMop.mass);
570: MatDiagonalScale(A,appctx->SEMop.mass,0);
571: VecReciprocal(appctx->SEMop.mass);
573: PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
574: return(0);
575: }
577: /* ------------------------------------------------------------------ */
578: /*
579: FormFunctionGradient - Evaluates the function and corresponding gradient.
581: Input Parameters:
582: tao - the Tao context
583: ic - the input vector
584: ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradientRoutine()
586: Output Parameters:
587: f - the newly evaluated function
588: G - the newly evaluated gradient
590: Notes:
592: The forward equation is
593: M u_t = F(U)
594: which is converted to
595: u_t = M^{-1} F(u)
596: in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is
597: M^{-1} J
598: where J is the Jacobian of F. Now the adjoint equation is
599: M v_t = J^T v
600: but TSAdjoint does not solve this since it can only solve the transposed system for the
601: Jacobian the user provided. Hence TSAdjoint solves
602: w_t = J^T M^{-1} w (where w = M v)
603: since there is no way to indicate the mass matrix as a separate entitity to TS. Thus one
604: must be careful in initializing the "adjoint equation" and using the result. This is
605: why
606: G = -2 M(u(T) - u_d)
607: below (instead of -2(u(T) - u_d)
610: */
611: PetscErrorCode FormFunctionGradient(Tao tao,Vec ic,PetscReal *f,Vec G,void *ctx)
612: {
613: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
614: PetscErrorCode ierr;
615: Vec temp;
618: TSSetTime(appctx->ts,0.0);
619: TSSetStepNumber(appctx->ts,0);
620: TSSetTimeStep(appctx->ts,appctx->initial_dt);
621: VecCopy(ic,appctx->dat.curr_sol);
623: TSSolve(appctx->ts,appctx->dat.curr_sol);
624: VecCopy(appctx->dat.curr_sol,appctx->dat.joe);
626: /* Compute the difference between the current ODE solution and target ODE solution */
627: VecWAXPY(G,-1.0,appctx->dat.curr_sol,appctx->dat.reference);
629: /* Compute the objective/cost function */
630: VecDuplicate(G,&temp);
631: VecPointwiseMult(temp,G,G);
632: VecDot(temp,appctx->SEMop.mass,f);
633: VecDestroy(&temp);
635: /* Compute initial conditions for the adjoint integration. See Notes above */
636: VecScale(G, -2.0);
637: VecPointwiseMult(G,G,appctx->SEMop.mass);
638: TSSetCostGradients(appctx->ts,1,&G,NULL);
640: TSAdjointSolve(appctx->ts);
641: /* VecPointwiseDivide(G,G,appctx->SEMop.mass);*/
642: return(0);
643: }
645: PetscErrorCode MonitorError(Tao tao,void *ctx)
646: {
647: AppCtx *appctx = (AppCtx*)ctx;
648: Vec temp,grad;
649: PetscReal nrm;
651: PetscInt its;
652: PetscReal fct,gnorm;
655: VecDuplicate(appctx->dat.ic,&temp);
656: VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution);
657: VecPointwiseMult(temp,temp,temp);
658: VecDot(temp,appctx->SEMop.mass,&nrm);
659: nrm = PetscSqrtReal(nrm);
660: TaoGetGradientVector(tao,&grad);
661: VecPointwiseMult(temp,temp,temp);
662: VecDot(temp,appctx->SEMop.mass,&gnorm);
663: gnorm = PetscSqrtReal(gnorm);
664: VecDestroy(&temp);
665: TaoGetIterationNumber(tao,&its);
666: TaoGetObjective(tao,&fct);
667: if (!its) {
668: PetscPrintf(PETSC_COMM_WORLD,"%% Iteration Error Objective Gradient-norm\n");
669: PetscPrintf(PETSC_COMM_WORLD,"history = [\n");
670: }
671: PetscPrintf(PETSC_COMM_WORLD,"%3D %g %g %g\n",its,(double)nrm,(double)fct,(double)gnorm);
672: return(0);
673: }
675: PetscErrorCode MonitorDestroy(void **ctx)
676: {
680: PetscPrintf(PETSC_COMM_WORLD,"];\n");
681: return(0);
682: }
685: /*TEST
687: build:
688: requires: !complex
690: test:
691: requires: !single
692: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
694: test:
695: suffix: cn
696: requires: !single
697: args: -ts_type cn -ts_dt .003 -pc_type lu -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
699: test:
700: suffix: 2
701: requires: !single
702: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -a .1 -tao_bqnls_mat_lmvm_scale_type none
705: TEST*/