Actual source code: ex50.c
petsc-3.11.4 2019-09-28
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/examples/tutorials/burgers_spectral.c
29: ------------------------------------------------------------------------- */
31: #include <petscts.h>
32: #include <petscgll.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; /* grid points per elements*/
43: PetscInt E; /* number of elements */
44: PetscReal tol_L2,tol_max; /* error norms */
45: PetscInt steps; /* number of timesteps */
46: PetscReal Tend; /* endtime */
47: PetscReal mu; /* viscosity */
48: PetscReal L; /* total length of domain */
49: PetscReal Le;
50: PetscReal Tadj;
51: } PetscParam;
53: typedef struct {
54: Vec grid; /* total grid */
55: Vec curr_sol;
56: } PetscData;
58: typedef struct {
59: Vec grid; /* total grid */
60: Vec mass; /* mass matrix for total integration */
61: Mat stiff; /* stifness matrix */
62: Mat keptstiff;
63: Mat grad;
64: PetscGLL gll;
65: } PetscSEMOperators;
67: typedef struct {
68: DM da; /* distributed array data structure */
69: PetscSEMOperators SEMop;
70: PetscParam param;
71: PetscData dat;
72: TS ts;
73: PetscReal initial_dt;
74: } AppCtx;
76: /*
77: User-defined routines
78: */
79: extern PetscErrorCode RHSMatrixLaplaciangllDM(TS,PetscReal,Vec,Mat,Mat,void*);
80: extern PetscErrorCode RHSMatrixAdvectiongllDM(TS,PetscReal,Vec,Mat,Mat,void*);
81: extern PetscErrorCode TrueSolution(TS,PetscReal,Vec,AppCtx*);
82: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
83: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
85: int main(int argc,char **argv)
86: {
87: AppCtx appctx; /* user-defined application context */
89: PetscInt i, xs, xm, ind, j, lenglob;
90: PetscReal x, *wrk_ptr1, *wrk_ptr2;
91: MatNullSpace nsp;
92: PetscMPIInt size;
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Initialize program and set problem parameters
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
101: /*initialize parameters */
102: appctx.param.N = 10; /* order of the spectral element */
103: appctx.param.E = 10; /* number of elements */
104: appctx.param.L = 4.0; /* length of the domain */
105: appctx.param.mu = 0.01; /* diffusion coefficient */
106: appctx.initial_dt = 5e-3;
107: appctx.param.steps = PETSC_MAX_INT;
108: appctx.param.Tend = 4;
110: PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL);
111: PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL);
112: PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL);
113: PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL);
114: appctx.param.Le = appctx.param.L/appctx.param.E;
116: MPI_Comm_size(PETSC_COMM_WORLD,&size);
117: if (appctx.param.E % size) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of elements must be divisible by number of processes");
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Create GLL data structures
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: PetscGLLCreate(appctx.param.N,PETSCGLL_VIA_LINEARALGEBRA,&appctx.SEMop.gll);
123: lenglob = appctx.param.E*(appctx.param.N-1);
125: /*
126: Create distributed array (DMDA) to manage parallel grid and vectors
127: and to set up the ghost point communication pattern. There are E*(Nl-1)+1
128: total grid values spread equally among all the processors, except first and last
129: */
131: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);
132: DMSetFromOptions(appctx.da);
133: DMSetUp(appctx.da);
134:
135: /*
136: Extract global and local vectors from DMDA; we use these to store the
137: approximate solution. Then duplicate these for remaining vectors that
138: have the same types.
139: */
141: DMCreateGlobalVector(appctx.da,&appctx.dat.curr_sol);
142: VecDuplicate(appctx.dat.curr_sol,&appctx.SEMop.grid);
143: VecDuplicate(appctx.dat.curr_sol,&appctx.SEMop.mass);
145: DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
146: DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
147: DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
148:
149: /* Compute function over the locally owned part of the grid */
150:
151: xs=xs/(appctx.param.N-1);
152: xm=xm/(appctx.param.N-1);
153:
154: /*
155: Build total grid and mass over entire mesh (multi-elemental)
156: */
158: for (i=xs; i<xs+xm; i++) {
159: for (j=0; j<appctx.param.N-1; j++) {
160: x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i;
161: ind=i*(appctx.param.N-1)+j;
162: wrk_ptr1[ind]=x;
163: wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
164: if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
165: }
166: }
167: DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
168: DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Create matrix data structure; set matrix evaluation routine.
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173: DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);
174: DMCreateMatrix(appctx.da,&appctx.SEMop.stiff);
175: DMCreateMatrix(appctx.da,&appctx.SEMop.grad);
176: /*
177: For linear problems with a time-dependent f(u,t) in the equation
178: u_t = f(u,t), the user provides the discretized right-hand-side
179: as a time-dependent matrix.
180: */
181: RHSMatrixLaplaciangllDM(appctx.ts,0.0,appctx.dat.curr_sol,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx);
182: RHSMatrixAdvectiongllDM(appctx.ts,0.0,appctx.dat.curr_sol,appctx.SEMop.grad,appctx.SEMop.grad,&appctx);
183: /*
184: For linear problems with a time-dependent f(u,t) in the equation
185: u_t = f(u,t), the user provides the discretized right-hand-side
186: as a time-dependent matrix.
187: */
188:
189: MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff);
191: /* attach the null space to the matrix, this probably is not needed but does no harm */
192: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
193: MatSetNullSpace(appctx.SEMop.stiff,nsp);
194: MatSetNullSpace(appctx.SEMop.keptstiff,nsp);
195: MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL);
196: MatNullSpaceDestroy(&nsp);
197: /* attach the null space to the matrix, this probably is not needed but does no harm */
198: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
199: MatSetNullSpace(appctx.SEMop.grad,nsp);
200: MatNullSpaceTest(nsp,appctx.SEMop.grad,NULL);
201: MatNullSpaceDestroy(&nsp);
203: /* Create the TS solver that solves the ODE and its adjoint; set its options */
204: TSCreate(PETSC_COMM_WORLD,&appctx.ts);
205: TSSetProblemType(appctx.ts,TS_NONLINEAR);
206: TSSetType(appctx.ts,TSRK);
207: TSSetDM(appctx.ts,appctx.da);
208: TSSetTime(appctx.ts,0.0);
209: TSSetTimeStep(appctx.ts,appctx.initial_dt);
210: TSSetMaxSteps(appctx.ts,appctx.param.steps);
211: TSSetMaxTime(appctx.ts,appctx.param.Tend);
212: TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
213: TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL);
214: TSSetSaveTrajectory(appctx.ts);
215: TSSetFromOptions(appctx.ts);
216: TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
217: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx);
219: /* Set Initial conditions for the problem */
220: TrueSolution(appctx.ts,0,appctx.dat.curr_sol,&appctx);
222: TSSetSolutionFunction(appctx.ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void *))TrueSolution,&appctx);
223: TSSetTime(appctx.ts,0.0);
224: TSSetStepNumber(appctx.ts,0);
226: TSSolve(appctx.ts,appctx.dat.curr_sol);
228: MatDestroy(&appctx.SEMop.stiff);
229: MatDestroy(&appctx.SEMop.keptstiff);
230: MatDestroy(&appctx.SEMop.grad);
231: VecDestroy(&appctx.SEMop.grid);
232: VecDestroy(&appctx.SEMop.mass);
233: VecDestroy(&appctx.dat.curr_sol);
234: PetscGLLDestroy(&appctx.SEMop.gll);
235: DMDestroy(&appctx.da);
236: TSDestroy(&appctx.ts);
238: /*
239: Always call PetscFinalize() before exiting a program. This routine
240: - finalizes the PETSc libraries as well as MPI
241: - provides summary and diagnostic information if certain runtime
242: options are chosen (e.g., -log_summary).
243: */
244: PetscFinalize();
245: return ierr;
246: }
248: /*
249: TrueSolution() computes the true solution for the PDE
251: Input Parameter:
252: u - uninitialized solution vector (global)
253: appctx - user-defined application context
255: Output Parameter:
256: u - vector with solution at initial time (global)
257: */
258: PetscErrorCode TrueSolution(TS ts, PetscReal t, Vec u,AppCtx *appctx)
259: {
260: PetscScalar *s;
261: const PetscScalar *xg;
262: PetscErrorCode ierr;
263: PetscInt i,xs,xn;
265: DMDAVecGetArray(appctx->da,u,&s);
266: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
267: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
268: for (i=xs; i<xs+xn; i++) {
269: 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));
270: }
271: DMDAVecRestoreArray(appctx->da,u,&s);
272: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
273: return 0;
274: }
276: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
277: {
279: AppCtx *appctx = (AppCtx*)ctx;
282: MatMult(appctx->SEMop.grad,globalin,globalout); /* grad u */
283: VecPointwiseMult(globalout,globalin,globalout); /* u grad u */
284: VecScale(globalout, -1.0);
285: MatMultAdd(appctx->SEMop.keptstiff,globalin,globalout,globalout);
286: return(0);
287: }
289: /*
291: K is the discretiziation of the Laplacian
292: G is the discretization of the gradient
294: Computes Jacobian of K u + diag(u) G u which is given by
295: K + diag(u)G + diag(Gu)
296: */
297: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx)
298: {
300: AppCtx *appctx = (AppCtx*)ctx;
301: Vec Gglobalin;
304: /* A = diag(u) G */
306: MatCopy(appctx->SEMop.grad,A,SAME_NONZERO_PATTERN);
307: MatDiagonalScale(A,globalin,NULL);
309: /* A = A + diag(Gu) */
310: VecDuplicate(globalin,&Gglobalin);
311: MatMult(appctx->SEMop.grad,globalin,Gglobalin);
312: MatDiagonalSet(A,Gglobalin,ADD_VALUES);
313: VecDestroy(&Gglobalin);
315: /* A = K - A */
316: MatScale(A,-1.0);
317: MatAXPY(A,0.0,appctx->SEMop.keptstiff,SAME_NONZERO_PATTERN);
318: return(0);
319: }
321: /* --------------------------------------------------------------------- */
323: #include "petscblaslapack.h"
324: /*
325: Matrix free operation of 1d Laplacian and Grad for GLL spectral elements
326: */
327: PetscErrorCode MatMult_Laplacian(Mat A,Vec x,Vec y)
328: {
329: AppCtx *appctx;
330: PetscErrorCode ierr;
331: PetscReal **temp,vv;
332: PetscInt i,j,xs,xn;
333: Vec xlocal,ylocal;
334: const PetscScalar *xl;
335: PetscScalar *yl;
336: PetscBLASInt _One = 1,n;
337: PetscScalar _DOne = 1;
339: MatShellGetContext(A,&appctx);
340: DMGetLocalVector(appctx->da,&xlocal);
341: DMGlobalToLocalBegin(appctx->da,x,INSERT_VALUES,xlocal);
342: DMGlobalToLocalEnd(appctx->da,x,INSERT_VALUES,xlocal);
343: DMGetLocalVector(appctx->da,&ylocal);
344: VecSet(ylocal,0.0);
345: PetscGLLElementLaplacianCreate(&appctx->SEMop.gll,&temp);
346: for (i=0; i<appctx->param.N; i++) {
347: vv =-appctx->param.mu*2.0/appctx->param.Le;
348: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
349: }
350: DMDAVecGetArrayRead(appctx->da,xlocal,(void*)&xl);
351: DMDAVecGetArray(appctx->da,ylocal,&yl);
352: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
353: PetscBLASIntCast(appctx->param.N,&n);
354: for (j=xs; j<xs+xn; j += appctx->param.N-1) {
355: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&_DOne,&temp[0][0],&n,&xl[j],&_One,&_DOne,&yl[j],&_One));
356: }
357: DMDAVecRestoreArrayRead(appctx->da,xlocal,(void*)&xl);
358: DMDAVecRestoreArray(appctx->da,ylocal,&yl);
359: PetscGLLElementLaplacianDestroy(&appctx->SEMop.gll,&temp);
360: VecSet(y,0.0);
361: DMLocalToGlobalBegin(appctx->da,ylocal,ADD_VALUES,y);
362: DMLocalToGlobalEnd(appctx->da,ylocal,ADD_VALUES,y);
363: DMRestoreLocalVector(appctx->da,&xlocal);
364: DMRestoreLocalVector(appctx->da,&ylocal);
365: VecPointwiseDivide(y,y,appctx->SEMop.mass);
366: return 0;
367: }
369: PetscErrorCode MatMult_Advection(Mat A,Vec x,Vec y)
370: {
371: AppCtx *appctx;
372: PetscErrorCode ierr;
373: PetscReal **temp;
374: PetscInt j,xs,xn;
375: Vec xlocal,ylocal;
376: const PetscScalar *xl;
377: PetscScalar *yl;
378: PetscBLASInt _One = 1,n;
379: PetscScalar _DOne = 1;
381: MatShellGetContext(A,&appctx);
382: DMGetLocalVector(appctx->da,&xlocal);
383: DMGlobalToLocalBegin(appctx->da,x,INSERT_VALUES,xlocal);
384: DMGlobalToLocalEnd(appctx->da,x,INSERT_VALUES,xlocal);
385: DMGetLocalVector(appctx->da,&ylocal);
386: VecSet(ylocal,0.0);
387: PetscGLLElementAdvectionCreate(&appctx->SEMop.gll,&temp);
388: DMDAVecGetArrayRead(appctx->da,xlocal,(void*)&xl);
389: DMDAVecGetArray(appctx->da,ylocal,&yl);
390: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
391: PetscBLASIntCast(appctx->param.N,&n);
392: for (j=xs; j<xs+xn; j += appctx->param.N-1) {
393: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&_DOne,&temp[0][0],&n,&xl[j],&_One,&_DOne,&yl[j],&_One));
394: }
395: DMDAVecRestoreArrayRead(appctx->da,xlocal,(void*)&xl);
396: DMDAVecRestoreArray(appctx->da,ylocal,&yl);
397: PetscGLLElementAdvectionDestroy(&appctx->SEMop.gll,&temp);
398: VecSet(y,0.0);
399: DMLocalToGlobalBegin(appctx->da,ylocal,ADD_VALUES,y);
400: DMLocalToGlobalEnd(appctx->da,ylocal,ADD_VALUES,y);
401: DMRestoreLocalVector(appctx->da,&xlocal);
402: DMRestoreLocalVector(appctx->da,&ylocal);
403: VecPointwiseDivide(y,y,appctx->SEMop.mass);
404: VecScale(y,-1.0);
405: return 0;
406: }
408: /*
409: RHSMatrixLaplacian - User-provided routine to compute the right-hand-side
410: matrix for the Laplacian operator
412: Input Parameters:
413: ts - the TS context
414: t - current time (ignored)
415: X - current solution (ignored)
416: dummy - optional user-defined context, as set by TSetRHSJacobian()
418: Output Parameters:
419: AA - Jacobian matrix
420: BB - optionally different matrix from which the preconditioner is built
421: str - flag indicating matrix structure
423: */
424: PetscErrorCode RHSMatrixLaplaciangllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
425: {
426: PetscReal **temp;
427: PetscReal vv;
428: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
430: PetscInt i,xs,xn,l,j;
431: PetscInt *rowsDM;
432: PetscBool flg = PETSC_FALSE;
434: PetscOptionsGetBool(NULL,NULL,"-gll_mf",&flg,NULL);
436: if (!flg) {
437: /*
438: Creates the element stiffness matrix for the given gll
439: */
440: PetscGLLElementLaplacianCreate(&appctx->SEMop.gll,&temp);
441: /* workarround for clang analyzer warning: Division by zero */
442: if (appctx->param.N <= 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Spectral element order should be > 1");
444: /* scale by the size of the element */
445: for (i=0; i<appctx->param.N; i++) {
446: vv=-appctx->param.mu*2.0/appctx->param.Le;
447: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
448: }
450: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
451: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
453: xs = xs/(appctx->param.N-1);
454: xn = xn/(appctx->param.N-1);
456: PetscMalloc1(appctx->param.N,&rowsDM);
457: /*
458: loop over local elements
459: */
460: for (j=xs; j<xs+xn; j++) {
461: for (l=0; l<appctx->param.N; l++) {
462: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
463: }
464: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
465: }
466: PetscFree(rowsDM);
467: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
468: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
469: VecReciprocal(appctx->SEMop.mass);
470: MatDiagonalScale(A,appctx->SEMop.mass,0);
471: VecReciprocal(appctx->SEMop.mass);
473: PetscGLLElementLaplacianDestroy(&appctx->SEMop.gll,&temp);
474: } else {
475: MatSetType(A,MATSHELL);
476: MatSetUp(A);
477: MatShellSetContext(A,appctx);
478: MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_Laplacian);
479: }
480: return 0;
481: }
483: /*
484: RHSMatrixAdvection - User-provided routine to compute the right-hand-side
485: matrix for the Advection (gradient) operator.
487: Input Parameters:
488: ts - the TS context
489: t - current time
490: global_in - global input vector
491: dummy - optional user-defined context, as set by TSetRHSJacobian()
493: Output Parameters:
494: AA - Jacobian matrix
495: BB - optionally different preconditioning matrix
496: str - flag indicating matrix structure
498: */
499: PetscErrorCode RHSMatrixAdvectiongllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
500: {
501: PetscReal **temp;
502: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
504: PetscInt xs,xn,l,j;
505: PetscInt *rowsDM;
506: PetscBool flg = PETSC_FALSE;
508: PetscOptionsGetBool(NULL,NULL,"-gll_mf",&flg,NULL);
510: if (!flg) {
511: /*
512: Creates the advection matrix for the given gll
513: */
514: PetscGLLElementAdvectionCreate(&appctx->SEMop.gll,&temp);
515: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
516: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
517: xs = xs/(appctx->param.N-1);
518: xn = xn/(appctx->param.N-1);
520: PetscMalloc1(appctx->param.N,&rowsDM);
521: for (j=xs; j<xs+xn; j++) {
522: for (l=0; l<appctx->param.N; l++) {
523: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
524: }
525: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
526: }
527: PetscFree(rowsDM);
528: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
529: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
531: VecReciprocal(appctx->SEMop.mass);
532: MatDiagonalScale(A,appctx->SEMop.mass,0);
533: VecReciprocal(appctx->SEMop.mass);
535: PetscGLLElementAdvectionDestroy(&appctx->SEMop.gll,&temp);
536: } else {
537: MatSetType(A,MATSHELL);
538: MatSetUp(A);
539: MatShellSetContext(A,appctx);
540: MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_Advection);
541: }
542: return 0;
543: }
545: /*TEST
547: build:
548: requires: !complex
550: test:
551: suffix: 1
552: requires: !single
554: test:
555: suffix: 2
556: nsize: 5
557: requires: !single
559: test:
560: suffix: 3
561: requires: !single
562: args: -ts_view -ts_type beuler -gll_mf -pc_type none -ts_max_steps 5 -ts_monitor_error
564: test:
565: suffix: 4
566: requires: !single
567: args: -ts_view -ts_type beuler -pc_type none -ts_max_steps 5 -ts_monitor_error
569: TEST*/