Actual source code: spectraladjointassimilation.c

petsc-master 2019-09-15
Report Typos and Errors

  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 seperate 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*/