Actual source code: ex50.c


  2: static char help[] ="Solves one dimensional Burger's equation compares with exact solution\n\n";

  4: /*

  6:     Not yet tested in parallel

  8: */
  9: /*
 10:    Concepts: TS^time-dependent nonlinear problems
 11:    Concepts: TS^Burger's equation
 12:    Processors: n
 13: */

 15: /* ------------------------------------------------------------------------

 17:    This program uses the one-dimensional Burger's equation
 18:        u_t = mu*u_xx - u u_x,
 19:    on the domain 0 <= x <= 1, with periodic boundary conditions

 21:    The operators are discretized with the spectral element method

 23:    See the paper PDE-CONSTRAINED OPTIMIZATION WITH SPECTRAL ELEMENTS USING PETSC AND TAO
 24:    by OANA MARIN, EMIL CONSTANTINESCU, AND BARRY SMITH for details on the exact solution
 25:    used

 27:    See src/tao/unconstrained/tutorials/burgers_spectral.c

 29:   ------------------------------------------------------------------------- */

 31: #include <petscts.h>
 32: #include <petscdt.h>
 33: #include <petscdraw.h>
 34: #include <petscdmda.h>

 36: /*
 37:    User-defined application context - contains data needed by the
 38:    application-provided call-back routines.
 39: */

 41: typedef struct {
 42:   PetscInt  n;                /* number of nodes */
 43:   PetscReal *nodes;           /* GLL nodes */
 44:   PetscReal *weights;         /* GLL weights */
 45: } PetscGLL;

 47: typedef struct {
 48:   PetscInt    N;             /* grid points per elements*/
 49:   PetscInt    E;              /* number of elements */
 50:   PetscReal   tol_L2,tol_max; /* error norms */
 51:   PetscInt    steps;          /* number of timesteps */
 52:   PetscReal   Tend;           /* endtime */
 53:   PetscReal   mu;             /* viscosity */
 54:   PetscReal   L;              /* total length of domain */
 55:   PetscReal   Le;
 56:   PetscReal   Tadj;
 57: } PetscParam;

 59: typedef struct {
 60:   Vec         grid;              /* total grid */
 61:   Vec         curr_sol;
 62: } PetscData;

 64: typedef struct {
 65:   Vec         grid;              /* total grid */
 66:   Vec         mass;              /* mass matrix for total integration */
 67:   Mat         stiff;             /* stifness matrix */
 68:   Mat         keptstiff;
 69:   Mat         grad;
 70:   PetscGLL    gll;
 71: } PetscSEMOperators;

 73: typedef struct {
 74:   DM                da;                /* distributed array data structure */
 75:   PetscSEMOperators SEMop;
 76:   PetscParam        param;
 77:   PetscData         dat;
 78:   TS                ts;
 79:   PetscReal         initial_dt;
 80: } AppCtx;

 82: /*
 83:    User-defined routines
 84: */
 85: extern PetscErrorCode RHSMatrixLaplaciangllDM(TS,PetscReal,Vec,Mat,Mat,void*);
 86: extern PetscErrorCode RHSMatrixAdvectiongllDM(TS,PetscReal,Vec,Mat,Mat,void*);
 87: extern PetscErrorCode TrueSolution(TS,PetscReal,Vec,AppCtx*);
 88: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
 89: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);

 91: int main(int argc,char **argv)
 92: {
 93:   AppCtx         appctx;                 /* user-defined application context */
 95:   PetscInt       i, xs, xm, ind, j, lenglob;
 96:   PetscReal      x, *wrk_ptr1, *wrk_ptr2;
 97:   MatNullSpace   nsp;
 98:   PetscMPIInt    size;

100:    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Initialize program and set problem parameters
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

105:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

107:   /*initialize parameters */
108:   appctx.param.N    = 10;  /* order of the spectral element */
109:   appctx.param.E    = 10;  /* number of elements */
110:   appctx.param.L    = 4.0;  /* length of the domain */
111:   appctx.param.mu   = 0.01; /* diffusion coefficient */
112:   appctx.initial_dt = 5e-3;
113:   appctx.param.steps = PETSC_MAX_INT;
114:   appctx.param.Tend  = 4;

116:   PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL);
117:   PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL);
118:   PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL);
119:   PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL);
120:   appctx.param.Le = appctx.param.L/appctx.param.E;

122:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
123:   if (appctx.param.E % size) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of elements must be divisible by number of processes");

125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:      Create GLL data structures
127:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128:   PetscMalloc2(appctx.param.N,&appctx.SEMop.gll.nodes,appctx.param.N,&appctx.SEMop.gll.weights);
129:   PetscDTGaussLobattoLegendreQuadrature(appctx.param.N,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
130:   appctx.SEMop.gll.n = appctx.param.N;
131:   lenglob  = appctx.param.E*(appctx.param.N-1);

133:   /*
134:      Create distributed array (DMDA) to manage parallel grid and vectors
135:      and to set up the ghost point communication pattern.  There are E*(Nl-1)+1
136:      total grid values spread equally among all the processors, except first and last
137:   */

139:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);
140:   DMSetFromOptions(appctx.da);
141:   DMSetUp(appctx.da);

143:   /*
144:      Extract global and local vectors from DMDA; we use these to store the
145:      approximate solution.  Then duplicate these for remaining vectors that
146:      have the same types.
147:   */

149:   DMCreateGlobalVector(appctx.da,&appctx.dat.curr_sol);
150:   VecDuplicate(appctx.dat.curr_sol,&appctx.SEMop.grid);
151:   VecDuplicate(appctx.dat.curr_sol,&appctx.SEMop.mass);

153:   DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
154:   DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
155:   DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);

157:   /* Compute function over the locally owned part of the grid */

159:     xs=xs/(appctx.param.N-1);
160:     xm=xm/(appctx.param.N-1);

162:   /*
163:      Build total grid and mass over entire mesh (multi-elemental)
164:   */

166:   for (i=xs; i<xs+xm; i++) {
167:     for (j=0; j<appctx.param.N-1; j++) {
168:       x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i;
169:       ind=i*(appctx.param.N-1)+j;
170:       wrk_ptr1[ind]=x;
171:       wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
172:       if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
173:     }
174:   }
175:   DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
176:   DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);

178:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179:    Create matrix data structure; set matrix evaluation routine.
180:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181:   DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);
182:   DMCreateMatrix(appctx.da,&appctx.SEMop.stiff);
183:   DMCreateMatrix(appctx.da,&appctx.SEMop.grad);
184:   /*
185:    For linear problems with a time-dependent f(u,t) in the equation
186:    u_t = f(u,t), the user provides the discretized right-hand-side
187:    as a time-dependent matrix.
188:    */
189:   RHSMatrixLaplaciangllDM(appctx.ts,0.0,appctx.dat.curr_sol,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx);
190:   RHSMatrixAdvectiongllDM(appctx.ts,0.0,appctx.dat.curr_sol,appctx.SEMop.grad,appctx.SEMop.grad,&appctx);
191:    /*
192:        For linear problems with a time-dependent f(u,t) in the equation
193:        u_t = f(u,t), the user provides the discretized right-hand-side
194:        as a time-dependent matrix.
195:     */

197:   MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff);

199:   /* attach the null space to the matrix, this probably is not needed but does no harm */
200:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
201:   MatSetNullSpace(appctx.SEMop.stiff,nsp);
202:   MatSetNullSpace(appctx.SEMop.keptstiff,nsp);
203:   MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL);
204:   MatNullSpaceDestroy(&nsp);
205:   /* attach the null space to the matrix, this probably is not needed but does no harm */
206:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
207:   MatSetNullSpace(appctx.SEMop.grad,nsp);
208:   MatNullSpaceTest(nsp,appctx.SEMop.grad,NULL);
209:   MatNullSpaceDestroy(&nsp);

211:   /* Create the TS solver that solves the ODE and its adjoint; set its options */
212:   TSCreate(PETSC_COMM_WORLD,&appctx.ts);
213:   TSSetProblemType(appctx.ts,TS_NONLINEAR);
214:   TSSetType(appctx.ts,TSRK);
215:   TSSetDM(appctx.ts,appctx.da);
216:   TSSetTime(appctx.ts,0.0);
217:   TSSetTimeStep(appctx.ts,appctx.initial_dt);
218:   TSSetMaxSteps(appctx.ts,appctx.param.steps);
219:   TSSetMaxTime(appctx.ts,appctx.param.Tend);
220:   TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
221:   TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL);
222:   TSSetSaveTrajectory(appctx.ts);
223:   TSSetFromOptions(appctx.ts);
224:   TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
225:   TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx);

227:   /* Set Initial conditions for the problem  */
228:   TrueSolution(appctx.ts,0,appctx.dat.curr_sol,&appctx);

230:   TSSetSolutionFunction(appctx.ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void *))TrueSolution,&appctx);
231:   TSSetTime(appctx.ts,0.0);
232:   TSSetStepNumber(appctx.ts,0);

234:   TSSolve(appctx.ts,appctx.dat.curr_sol);

236:   MatDestroy(&appctx.SEMop.stiff);
237:   MatDestroy(&appctx.SEMop.keptstiff);
238:   MatDestroy(&appctx.SEMop.grad);
239:   VecDestroy(&appctx.SEMop.grid);
240:   VecDestroy(&appctx.SEMop.mass);
241:   VecDestroy(&appctx.dat.curr_sol);
242:   PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights);
243:   DMDestroy(&appctx.da);
244:   TSDestroy(&appctx.ts);

246:   /*
247:      Always call PetscFinalize() before exiting a program.  This routine
248:        - finalizes the PETSc libraries as well as MPI
249:        - provides summary and diagnostic information if certain runtime
250:          options are chosen (e.g., -log_summary).
251:   */
252:     PetscFinalize();
253:     return ierr;
254: }

256: /*
257:    TrueSolution() computes the true solution for the PDE

259:    Input Parameter:
260:    u - uninitialized solution vector (global)
261:    appctx - user-defined application context

263:    Output Parameter:
264:    u - vector with solution at initial time (global)
265: */
266: PetscErrorCode TrueSolution(TS ts, PetscReal t, Vec u,AppCtx *appctx)
267: {
268:   PetscScalar       *s;
269:   const PetscScalar *xg;
270:   PetscErrorCode    ierr;
271:   PetscInt          i,xs,xn;

273:   DMDAVecGetArray(appctx->da,u,&s);
274:   DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
275:   DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
276:   for (i=xs; i<xs+xn; i++) {
277:     s[i]=2.0*appctx->param.mu*PETSC_PI*PetscSinScalar(PETSC_PI*xg[i])*PetscExpReal(-appctx->param.mu*PETSC_PI*PETSC_PI*t)/(2.0+PetscCosScalar(PETSC_PI*xg[i])*PetscExpReal(-appctx->param.mu*PETSC_PI*PETSC_PI*t));
278:   }
279:   DMDAVecRestoreArray(appctx->da,u,&s);
280:   DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
281:   return 0;
282: }

284: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
285: {
287:   AppCtx          *appctx = (AppCtx*)ctx;

290:   MatMult(appctx->SEMop.grad,globalin,globalout); /* grad u */
291:   VecPointwiseMult(globalout,globalin,globalout); /* u grad u */
292:   VecScale(globalout, -1.0);
293:   MatMultAdd(appctx->SEMop.keptstiff,globalin,globalout,globalout);
294:   return(0);
295: }

297: /*

299:       K is the discretiziation of the Laplacian
300:       G is the discretization of the gradient

302:       Computes Jacobian of      K u + diag(u) G u   which is given by
303:               K   + diag(u)G + diag(Gu)
304: */
305: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx)
306: {
308:   AppCtx         *appctx = (AppCtx*)ctx;
309:   Vec            Gglobalin;

312:   /*    A = diag(u) G */

314:   MatCopy(appctx->SEMop.grad,A,SAME_NONZERO_PATTERN);
315:   MatDiagonalScale(A,globalin,NULL);

317:   /*    A  = A + diag(Gu) */
318:   VecDuplicate(globalin,&Gglobalin);
319:   MatMult(appctx->SEMop.grad,globalin,Gglobalin);
320:   MatDiagonalSet(A,Gglobalin,ADD_VALUES);
321:   VecDestroy(&Gglobalin);

323:   /*   A  = K - A    */
324:   MatScale(A,-1.0);
325:   MatAXPY(A,0.0,appctx->SEMop.keptstiff,SAME_NONZERO_PATTERN);
326:   return(0);
327: }

329: /* --------------------------------------------------------------------- */

331: #include "petscblaslapack.h"
332: /*
333:      Matrix free operation of 1d Laplacian and Grad for GLL spectral elements
334: */
335: PetscErrorCode MatMult_Laplacian(Mat A,Vec x,Vec y)
336: {
337:   AppCtx            *appctx;
338:   PetscErrorCode    ierr;
339:   PetscReal         **temp,vv;
340:   PetscInt          i,j,xs,xn;
341:   Vec               xlocal,ylocal;
342:   const PetscScalar *xl;
343:   PetscScalar       *yl;
344:   PetscBLASInt      _One = 1,n;
345:   PetscScalar       _DOne = 1;

347:   MatShellGetContext(A,&appctx);
348:   DMGetLocalVector(appctx->da,&xlocal);
349:   DMGlobalToLocalBegin(appctx->da,x,INSERT_VALUES,xlocal);
350:   DMGlobalToLocalEnd(appctx->da,x,INSERT_VALUES,xlocal);
351:   DMGetLocalVector(appctx->da,&ylocal);
352:   VecSet(ylocal,0.0);
353:   PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
354:   for (i=0; i<appctx->param.N; i++) {
355:     vv =-appctx->param.mu*2.0/appctx->param.Le;
356:     for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
357:   }
358:   DMDAVecGetArrayRead(appctx->da,xlocal,(void*)&xl);
359:   DMDAVecGetArray(appctx->da,ylocal,&yl);
360:   DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
361:   PetscBLASIntCast(appctx->param.N,&n);
362:   for (j=xs; j<xs+xn; j += appctx->param.N-1) {
363:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&_DOne,&temp[0][0],&n,&xl[j],&_One,&_DOne,&yl[j],&_One));
364:   }
365:   DMDAVecRestoreArrayRead(appctx->da,xlocal,(void*)&xl);
366:   DMDAVecRestoreArray(appctx->da,ylocal,&yl);
367:   PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
368:   VecSet(y,0.0);
369:   DMLocalToGlobalBegin(appctx->da,ylocal,ADD_VALUES,y);
370:   DMLocalToGlobalEnd(appctx->da,ylocal,ADD_VALUES,y);
371:   DMRestoreLocalVector(appctx->da,&xlocal);
372:   DMRestoreLocalVector(appctx->da,&ylocal);
373:   VecPointwiseDivide(y,y,appctx->SEMop.mass);
374:   return 0;
375: }

377: PetscErrorCode MatMult_Advection(Mat A,Vec x,Vec y)
378: {
379:   AppCtx            *appctx;
380:   PetscErrorCode    ierr;
381:   PetscReal         **temp;
382:   PetscInt          j,xs,xn;
383:   Vec               xlocal,ylocal;
384:   const PetscScalar *xl;
385:   PetscScalar       *yl;
386:   PetscBLASInt      _One = 1,n;
387:   PetscScalar       _DOne = 1;

389:   MatShellGetContext(A,&appctx);
390:   DMGetLocalVector(appctx->da,&xlocal);
391:   DMGlobalToLocalBegin(appctx->da,x,INSERT_VALUES,xlocal);
392:   DMGlobalToLocalEnd(appctx->da,x,INSERT_VALUES,xlocal);
393:   DMGetLocalVector(appctx->da,&ylocal);
394:   VecSet(ylocal,0.0);
395:   PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
396:   DMDAVecGetArrayRead(appctx->da,xlocal,(void*)&xl);
397:   DMDAVecGetArray(appctx->da,ylocal,&yl);
398:   DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
399:   PetscBLASIntCast(appctx->param.N,&n);
400:   for (j=xs; j<xs+xn; j += appctx->param.N-1) {
401:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&_DOne,&temp[0][0],&n,&xl[j],&_One,&_DOne,&yl[j],&_One));
402:   }
403:   DMDAVecRestoreArrayRead(appctx->da,xlocal,(void*)&xl);
404:   DMDAVecRestoreArray(appctx->da,ylocal,&yl);
405:   PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
406:   VecSet(y,0.0);
407:   DMLocalToGlobalBegin(appctx->da,ylocal,ADD_VALUES,y);
408:   DMLocalToGlobalEnd(appctx->da,ylocal,ADD_VALUES,y);
409:   DMRestoreLocalVector(appctx->da,&xlocal);
410:   DMRestoreLocalVector(appctx->da,&ylocal);
411:   VecPointwiseDivide(y,y,appctx->SEMop.mass);
412:   VecScale(y,-1.0);
413:   return 0;
414: }

416: /*
417:    RHSMatrixLaplacian - User-provided routine to compute the right-hand-side
418:    matrix for the Laplacian operator

420:    Input Parameters:
421:    ts - the TS context
422:    t - current time  (ignored)
423:    X - current solution (ignored)
424:    dummy - optional user-defined context, as set by TSetRHSJacobian()

426:    Output Parameters:
427:    AA - Jacobian matrix
428:    BB - optionally different matrix from which the preconditioner is built
429:    str - flag indicating matrix structure

431: */
432: PetscErrorCode RHSMatrixLaplaciangllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
433: {
434:   PetscReal      **temp;
435:   PetscReal      vv;
436:   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */
438:   PetscInt       i,xs,xn,l,j;
439:   PetscInt       *rowsDM;
440:   PetscBool      flg = PETSC_FALSE;

442:   PetscOptionsGetBool(NULL,NULL,"-gll_mf",&flg,NULL);

444:   if (!flg) {
445:     /*
446:      Creates the element stiffness matrix for the given gll
447:      */
448:     PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
449:     /* workarround for clang analyzer warning: Division by zero */
450:     if (appctx->param.N <= 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Spectral element order should be > 1");

452:     /* scale by the size of the element */
453:     for (i=0; i<appctx->param.N; i++) {
454:       vv=-appctx->param.mu*2.0/appctx->param.Le;
455:       for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
456:     }

458:     MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
459:     DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);

461:     xs   = xs/(appctx->param.N-1);
462:     xn   = xn/(appctx->param.N-1);

464:     PetscMalloc1(appctx->param.N,&rowsDM);
465:     /*
466:      loop over local elements
467:      */
468:     for (j=xs; j<xs+xn; j++) {
469:       for (l=0; l<appctx->param.N; l++) {
470:         rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
471:       }
472:       MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
473:     }
474:     PetscFree(rowsDM);
475:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
476:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
477:     VecReciprocal(appctx->SEMop.mass);
478:     MatDiagonalScale(A,appctx->SEMop.mass,0);
479:     VecReciprocal(appctx->SEMop.mass);

481:     PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
482:   } else {
483:     MatSetType(A,MATSHELL);
484:     MatSetUp(A);
485:     MatShellSetContext(A,appctx);
486:     MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_Laplacian);
487:   }
488:   return 0;
489: }

491: /*
492:    RHSMatrixAdvection - User-provided routine to compute the right-hand-side
493:    matrix for the Advection (gradient) operator.

495:    Input Parameters:
496:    ts - the TS context
497:    t - current time
498:    global_in - global input vector
499:    dummy - optional user-defined context, as set by TSetRHSJacobian()

501:    Output Parameters:
502:    AA - Jacobian matrix
503:    BB - optionally different preconditioning matrix
504:    str - flag indicating matrix structure

506: */
507: PetscErrorCode RHSMatrixAdvectiongllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
508: {
509:   PetscReal      **temp;
510:   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */
512:   PetscInt       xs,xn,l,j;
513:   PetscInt       *rowsDM;
514:   PetscBool      flg = PETSC_FALSE;

516:   PetscOptionsGetBool(NULL,NULL,"-gll_mf",&flg,NULL);

518:   if (!flg) {
519:     /*
520:      Creates the advection matrix for the given gll
521:      */
522:     PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
523:     MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
524:     DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
525:     xs   = xs/(appctx->param.N-1);
526:     xn   = xn/(appctx->param.N-1);

528:     PetscMalloc1(appctx->param.N,&rowsDM);
529:     for (j=xs; j<xs+xn; j++) {
530:       for (l=0; l<appctx->param.N; l++) {
531:         rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
532:       }
533:       MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
534:     }
535:     PetscFree(rowsDM);
536:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
537:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

539:     VecReciprocal(appctx->SEMop.mass);
540:     MatDiagonalScale(A,appctx->SEMop.mass,0);
541:     VecReciprocal(appctx->SEMop.mass);

543:     PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp);
544:   } else {
545:     MatSetType(A,MATSHELL);
546:     MatSetUp(A);
547:     MatShellSetContext(A,appctx);
548:     MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_Advection);
549:   }
550:   return 0;
551: }

553: /*TEST

555:     build:
556:       requires: !complex

558:     test:
559:       suffix: 1
560:       requires: !single

562:     test:
563:       suffix: 2
564:       nsize: 5
565:       requires: !single

567:     test:
568:       suffix: 3
569:       requires: !single
570:       args: -ts_view  -ts_type beuler -gll_mf -pc_type none -ts_max_steps 5 -ts_monitor_error

572:     test:
573:       suffix: 4
574:       requires: !single
575:       args: -ts_view  -ts_type beuler  -pc_type none -ts_max_steps 5 -ts_monitor_error

577: TEST*/