Actual source code: ex50.c
1: static char help[] = "Solves one dimensional Burger's equation compares with exact solution\n\n";
3: /*
4: Not yet tested in parallel
5: */
7: /* ------------------------------------------------------------------------
9: This program uses the one-dimensional Burger's equation
10: u_t = mu*u_xx - u u_x,
11: on the domain 0 <= x <= 1, with periodic boundary conditions
13: The operators are discretized with the spectral element method
15: See the paper PDE-CONSTRAINED OPTIMIZATION WITH SPECTRAL ELEMENTS USING PETSC AND TAO
16: by OANA MARIN, EMIL CONSTANTINESCU, AND BARRY SMITH for details on the exact solution
17: used
19: See src/tao/unconstrained/tutorials/burgers_spectral.c
21: ------------------------------------------------------------------------- */
23: #include <petscts.h>
24: #include <petscdt.h>
25: #include <petscdmda.h>
27: /*
28: User-defined application context - contains data needed by the
29: application-provided call-back routines.
30: */
32: typedef struct {
33: PetscInt n; /* number of nodes */
34: PetscReal *nodes; /* GLL nodes */
35: PetscReal *weights; /* GLL weights */
36: } PetscGLL;
38: typedef struct {
39: PetscInt N; /* grid points per elements*/
40: PetscInt E; /* number of elements */
41: PetscReal tol_L2, tol_max; /* error norms */
42: PetscInt steps; /* number of timesteps */
43: PetscReal Tend; /* endtime */
44: PetscReal mu; /* viscosity */
45: PetscReal L; /* total length of domain */
46: PetscReal Le;
47: PetscReal Tadj;
48: } PetscParam;
50: typedef struct {
51: Vec grid; /* total grid */
52: Vec curr_sol;
53: } PetscData;
55: typedef struct {
56: Vec grid; /* total grid */
57: Vec mass; /* mass matrix for total integration */
58: Mat stiff; /* stiffness matrix */
59: Mat keptstiff;
60: Mat grad;
61: PetscGLL gll;
62: } PetscSEMOperators;
64: typedef struct {
65: DM da; /* distributed array data structure */
66: PetscSEMOperators SEMop;
67: PetscParam param;
68: PetscData dat;
69: TS ts;
70: PetscReal initial_dt;
71: } AppCtx;
73: /*
74: User-defined routines
75: */
76: extern PetscErrorCode RHSMatrixLaplaciangllDM(TS, PetscReal, Vec, Mat, Mat, void *);
77: extern PetscErrorCode RHSMatrixAdvectiongllDM(TS, PetscReal, Vec, Mat, Mat, void *);
78: extern PetscErrorCode TrueSolution(TS, PetscReal, Vec, AppCtx *);
79: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
80: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
82: int main(int argc, char **argv)
83: {
84: AppCtx appctx; /* user-defined application context */
85: PetscInt i, xs, xm, ind, j, lenglob;
86: PetscReal x, *wrk_ptr1, *wrk_ptr2;
87: MatNullSpace nsp;
88: PetscMPIInt size;
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Initialize program and set problem parameters
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: PetscFunctionBeginUser;
94: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
96: /*initialize parameters */
97: appctx.param.N = 10; /* order of the spectral element */
98: appctx.param.E = 10; /* number of elements */
99: appctx.param.L = 4.0; /* length of the domain */
100: appctx.param.mu = 0.01; /* diffusion coefficient */
101: appctx.initial_dt = 5e-3;
102: appctx.param.steps = PETSC_INT_MAX;
103: appctx.param.Tend = 4;
105: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &appctx.param.N, NULL));
106: PetscCall(PetscOptionsGetInt(NULL, NULL, "-E", &appctx.param.E, NULL));
107: PetscCall(PetscOptionsGetReal(NULL, NULL, "-Tend", &appctx.param.Tend, NULL));
108: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &appctx.param.mu, NULL));
109: appctx.param.Le = appctx.param.L / appctx.param.E;
111: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
112: PetscCheck((appctx.param.E % size) == 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of elements must be divisible by number of processes");
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Create GLL data structures
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: PetscCall(PetscMalloc2(appctx.param.N, &appctx.SEMop.gll.nodes, appctx.param.N, &appctx.SEMop.gll.weights));
118: PetscCall(PetscDTGaussLobattoLegendreQuadrature(appctx.param.N, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, appctx.SEMop.gll.nodes, appctx.SEMop.gll.weights));
119: appctx.SEMop.gll.n = appctx.param.N;
120: lenglob = appctx.param.E * (appctx.param.N - 1);
122: /*
123: Create distributed array (DMDA) to manage parallel grid and vectors
124: and to set up the ghost point communication pattern. There are E*(Nl-1)+1
125: total grid values spread equally among all the processors, except first and last
126: */
128: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, lenglob, 1, 1, NULL, &appctx.da));
129: PetscCall(DMSetFromOptions(appctx.da));
130: PetscCall(DMSetUp(appctx.da));
132: /*
133: Extract global and local vectors from DMDA; we use these to store the
134: approximate solution. Then duplicate these for remaining vectors that
135: have the same types.
136: */
138: PetscCall(DMCreateGlobalVector(appctx.da, &appctx.dat.curr_sol));
139: PetscCall(VecDuplicate(appctx.dat.curr_sol, &appctx.SEMop.grid));
140: PetscCall(VecDuplicate(appctx.dat.curr_sol, &appctx.SEMop.mass));
142: PetscCall(DMDAGetCorners(appctx.da, &xs, NULL, NULL, &xm, NULL, NULL));
143: PetscCall(DMDAVecGetArray(appctx.da, appctx.SEMop.grid, &wrk_ptr1));
144: PetscCall(DMDAVecGetArray(appctx.da, appctx.SEMop.mass, &wrk_ptr2));
146: /* Compute function over the locally owned part of the grid */
148: xs = xs / (appctx.param.N - 1);
149: xm = xm / (appctx.param.N - 1);
151: /*
152: Build total grid and mass over entire mesh (multi-elemental)
153: */
155: for (i = xs; i < xs + xm; i++) {
156: for (j = 0; j < appctx.param.N - 1; j++) {
157: x = (appctx.param.Le / 2.0) * (appctx.SEMop.gll.nodes[j] + 1.0) + appctx.param.Le * i;
158: ind = i * (appctx.param.N - 1) + j;
159: wrk_ptr1[ind] = x;
160: wrk_ptr2[ind] = .5 * appctx.param.Le * appctx.SEMop.gll.weights[j];
161: if (j == 0) wrk_ptr2[ind] += .5 * appctx.param.Le * appctx.SEMop.gll.weights[j];
162: }
163: }
164: PetscCall(DMDAVecRestoreArray(appctx.da, appctx.SEMop.grid, &wrk_ptr1));
165: PetscCall(DMDAVecRestoreArray(appctx.da, appctx.SEMop.mass, &wrk_ptr2));
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Create matrix data structure; set matrix evaluation routine.
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: PetscCall(DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE));
171: PetscCall(DMCreateMatrix(appctx.da, &appctx.SEMop.stiff));
172: PetscCall(DMCreateMatrix(appctx.da, &appctx.SEMop.grad));
173: /*
174: For linear problems with a time-dependent f(u,t) in the equation
175: u_t = f(u,t), the user provides the discretized right-hand side
176: as a time-dependent matrix.
177: */
178: PetscCall(RHSMatrixLaplaciangllDM(appctx.ts, 0.0, appctx.dat.curr_sol, appctx.SEMop.stiff, appctx.SEMop.stiff, &appctx));
179: PetscCall(RHSMatrixAdvectiongllDM(appctx.ts, 0.0, appctx.dat.curr_sol, appctx.SEMop.grad, appctx.SEMop.grad, &appctx));
180: /*
181: For linear problems with a time-dependent f(u,t) in the equation
182: u_t = f(u,t), the user provides the discretized right-hand side
183: as a time-dependent matrix.
184: */
186: PetscCall(MatDuplicate(appctx.SEMop.stiff, MAT_COPY_VALUES, &appctx.SEMop.keptstiff));
188: /* attach the null space to the matrix, this probably is not needed but does no harm */
189: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp));
190: PetscCall(MatSetNullSpace(appctx.SEMop.stiff, nsp));
191: PetscCall(MatSetNullSpace(appctx.SEMop.keptstiff, nsp));
192: PetscCall(MatNullSpaceTest(nsp, appctx.SEMop.stiff, NULL));
193: PetscCall(MatNullSpaceDestroy(&nsp));
194: /* attach the null space to the matrix, this probably is not needed but does no harm */
195: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp));
196: PetscCall(MatSetNullSpace(appctx.SEMop.grad, nsp));
197: PetscCall(MatNullSpaceTest(nsp, appctx.SEMop.grad, NULL));
198: PetscCall(MatNullSpaceDestroy(&nsp));
200: /* Create the TS solver that solves the ODE and its adjoint; set its options */
201: PetscCall(TSCreate(PETSC_COMM_WORLD, &appctx.ts));
202: PetscCall(TSSetProblemType(appctx.ts, TS_NONLINEAR));
203: PetscCall(TSSetType(appctx.ts, TSRK));
204: PetscCall(TSSetDM(appctx.ts, appctx.da));
205: PetscCall(TSSetTime(appctx.ts, 0.0));
206: PetscCall(TSSetTimeStep(appctx.ts, appctx.initial_dt));
207: PetscCall(TSSetMaxSteps(appctx.ts, appctx.param.steps));
208: PetscCall(TSSetMaxTime(appctx.ts, appctx.param.Tend));
209: PetscCall(TSSetExactFinalTime(appctx.ts, TS_EXACTFINALTIME_MATCHSTEP));
210: PetscCall(TSSetTolerances(appctx.ts, 1e-7, NULL, 1e-7, NULL));
211: PetscCall(TSSetSaveTrajectory(appctx.ts));
212: PetscCall(TSSetFromOptions(appctx.ts));
213: PetscCall(TSSetRHSFunction(appctx.ts, NULL, RHSFunction, &appctx));
214: PetscCall(TSSetRHSJacobian(appctx.ts, appctx.SEMop.stiff, appctx.SEMop.stiff, RHSJacobian, &appctx));
216: /* Set Initial conditions for the problem */
217: PetscCall(TrueSolution(appctx.ts, 0, appctx.dat.curr_sol, &appctx));
219: PetscCall(TSSetSolutionFunction(appctx.ts, (PetscErrorCode (*)(TS, PetscReal, Vec, void *))TrueSolution, &appctx));
220: PetscCall(TSSetTime(appctx.ts, 0.0));
221: PetscCall(TSSetStepNumber(appctx.ts, 0));
223: PetscCall(TSSolve(appctx.ts, appctx.dat.curr_sol));
225: PetscCall(MatDestroy(&appctx.SEMop.stiff));
226: PetscCall(MatDestroy(&appctx.SEMop.keptstiff));
227: PetscCall(MatDestroy(&appctx.SEMop.grad));
228: PetscCall(VecDestroy(&appctx.SEMop.grid));
229: PetscCall(VecDestroy(&appctx.SEMop.mass));
230: PetscCall(VecDestroy(&appctx.dat.curr_sol));
231: PetscCall(PetscFree2(appctx.SEMop.gll.nodes, appctx.SEMop.gll.weights));
232: PetscCall(DMDestroy(&appctx.da));
233: PetscCall(TSDestroy(&appctx.ts));
235: /*
236: Always call PetscFinalize() before exiting a program. This routine
237: - finalizes the PETSc libraries as well as MPI
238: - provides summary and diagnostic information if certain runtime
239: options are chosen (e.g., -log_view).
240: */
241: PetscCall(PetscFinalize());
242: return 0;
243: }
245: /*
246: TrueSolution() computes the true solution for the PDE
248: Input Parameter:
249: u - uninitialized solution vector (global)
250: appctx - user-defined application context
252: Output Parameter:
253: u - vector with solution at initial time (global)
254: */
255: PetscErrorCode TrueSolution(TS ts, PetscReal t, Vec u, AppCtx *appctx)
256: {
257: PetscScalar *s;
258: const PetscScalar *xg;
259: PetscInt i, xs, xn;
261: PetscFunctionBeginUser;
262: PetscCall(DMDAVecGetArray(appctx->da, u, &s));
263: PetscCall(DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
264: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
265: for (i = xs; i < xs + xn; i++) {
266: 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));
267: }
268: PetscCall(DMDAVecRestoreArray(appctx->da, u, &s));
269: PetscCall(DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, PetscCtx ctx)
274: {
275: AppCtx *appctx = (AppCtx *)ctx;
277: PetscFunctionBeginUser;
278: PetscCall(MatMult(appctx->SEMop.grad, globalin, globalout)); /* grad u */
279: PetscCall(VecPointwiseMult(globalout, globalin, globalout)); /* u grad u */
280: PetscCall(VecScale(globalout, -1.0));
281: PetscCall(MatMultAdd(appctx->SEMop.keptstiff, globalin, globalout, globalout));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*
287: K is the discretiziation of the Laplacian
288: G is the discretization of the gradient
290: Computes Jacobian of K u + diag(u) G u which is given by
291: K + diag(u)G + diag(Gu)
292: */
293: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec globalin, Mat A, Mat B, PetscCtx ctx)
294: {
295: AppCtx *appctx = (AppCtx *)ctx;
296: Vec Gglobalin;
298: PetscFunctionBeginUser;
299: /* A = diag(u) G */
301: PetscCall(MatCopy(appctx->SEMop.grad, A, SAME_NONZERO_PATTERN));
302: PetscCall(MatDiagonalScale(A, globalin, NULL));
304: /* A = A + diag(Gu) */
305: PetscCall(VecDuplicate(globalin, &Gglobalin));
306: PetscCall(MatMult(appctx->SEMop.grad, globalin, Gglobalin));
307: PetscCall(MatDiagonalSet(A, Gglobalin, ADD_VALUES));
308: PetscCall(VecDestroy(&Gglobalin));
310: /* A = K - A */
311: PetscCall(MatScale(A, -1.0));
312: PetscCall(MatAXPY(A, 0.0, appctx->SEMop.keptstiff, SAME_NONZERO_PATTERN));
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: #include <petscblaslapack.h>
317: /*
318: Matrix free operation of 1d Laplacian and Grad for GLL spectral elements
319: */
320: PetscErrorCode MatMult_Laplacian(Mat A, Vec x, Vec y)
321: {
322: AppCtx *appctx;
323: PetscReal **temp, vv;
324: PetscInt i, j, xs, xn;
325: Vec xlocal, ylocal;
326: const PetscScalar *xl;
327: PetscScalar *yl;
328: PetscBLASInt _One = 1, n;
329: PetscScalar _DOne = 1;
331: PetscFunctionBeginUser;
332: PetscCall(MatShellGetContext(A, &appctx));
333: PetscCall(DMGetLocalVector(appctx->da, &xlocal));
334: PetscCall(DMGlobalToLocalBegin(appctx->da, x, INSERT_VALUES, xlocal));
335: PetscCall(DMGlobalToLocalEnd(appctx->da, x, INSERT_VALUES, xlocal));
336: PetscCall(DMGetLocalVector(appctx->da, &ylocal));
337: PetscCall(VecSet(ylocal, 0.0));
338: PetscCall(PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
339: for (i = 0; i < appctx->param.N; i++) {
340: vv = -appctx->param.mu * 2.0 / appctx->param.Le;
341: for (j = 0; j < appctx->param.N; j++) temp[i][j] = temp[i][j] * vv;
342: }
343: PetscCall(DMDAVecGetArrayRead(appctx->da, xlocal, (void *)&xl));
344: PetscCall(DMDAVecGetArray(appctx->da, ylocal, &yl));
345: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
346: PetscCall(PetscBLASIntCast(appctx->param.N, &n));
347: for (j = xs; j < xs + xn; j += appctx->param.N - 1) PetscCallBLAS("BLASgemv", BLASgemv_("N", &n, &n, &_DOne, &temp[0][0], &n, &xl[j], &_One, &_DOne, &yl[j], &_One));
348: PetscCall(DMDAVecRestoreArrayRead(appctx->da, xlocal, (void *)&xl));
349: PetscCall(DMDAVecRestoreArray(appctx->da, ylocal, &yl));
350: PetscCall(PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
351: PetscCall(VecSet(y, 0.0));
352: PetscCall(DMLocalToGlobalBegin(appctx->da, ylocal, ADD_VALUES, y));
353: PetscCall(DMLocalToGlobalEnd(appctx->da, ylocal, ADD_VALUES, y));
354: PetscCall(DMRestoreLocalVector(appctx->da, &xlocal));
355: PetscCall(DMRestoreLocalVector(appctx->da, &ylocal));
356: PetscCall(VecPointwiseDivide(y, y, appctx->SEMop.mass));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: PetscErrorCode MatMult_Advection(Mat A, Vec x, Vec y)
361: {
362: AppCtx *appctx;
363: PetscReal **temp;
364: PetscInt j, xs, xn;
365: Vec xlocal, ylocal;
366: const PetscScalar *xl;
367: PetscScalar *yl;
368: PetscBLASInt _One = 1, n;
369: PetscScalar _DOne = 1;
371: PetscFunctionBeginUser;
372: PetscCall(MatShellGetContext(A, &appctx));
373: PetscCall(DMGetLocalVector(appctx->da, &xlocal));
374: PetscCall(DMGlobalToLocalBegin(appctx->da, x, INSERT_VALUES, xlocal));
375: PetscCall(DMGlobalToLocalEnd(appctx->da, x, INSERT_VALUES, xlocal));
376: PetscCall(DMGetLocalVector(appctx->da, &ylocal));
377: PetscCall(VecSet(ylocal, 0.0));
378: PetscCall(PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
379: PetscCall(DMDAVecGetArrayRead(appctx->da, xlocal, (void *)&xl));
380: PetscCall(DMDAVecGetArray(appctx->da, ylocal, &yl));
381: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
382: PetscCall(PetscBLASIntCast(appctx->param.N, &n));
383: for (j = xs; j < xs + xn; j += appctx->param.N - 1) PetscCallBLAS("BLASgemv", BLASgemv_("N", &n, &n, &_DOne, &temp[0][0], &n, &xl[j], &_One, &_DOne, &yl[j], &_One));
384: PetscCall(DMDAVecRestoreArrayRead(appctx->da, xlocal, (void *)&xl));
385: PetscCall(DMDAVecRestoreArray(appctx->da, ylocal, &yl));
386: PetscCall(PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
387: PetscCall(VecSet(y, 0.0));
388: PetscCall(DMLocalToGlobalBegin(appctx->da, ylocal, ADD_VALUES, y));
389: PetscCall(DMLocalToGlobalEnd(appctx->da, ylocal, ADD_VALUES, y));
390: PetscCall(DMRestoreLocalVector(appctx->da, &xlocal));
391: PetscCall(DMRestoreLocalVector(appctx->da, &ylocal));
392: PetscCall(VecPointwiseDivide(y, y, appctx->SEMop.mass));
393: PetscCall(VecScale(y, -1.0));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: /*
398: RHSMatrixLaplacian - User-provided routine to compute the right-hand-side
399: matrix for the Laplacian operator
401: Input Parameters:
402: ts - the TS context
403: t - current time (ignored)
404: X - current solution (ignored)
405: dummy - optional user-defined context, as set by TSetRHSJacobian()
407: Output Parameters:
408: AA - Jacobian matrix
409: BB - optionally different matrix from which the preconditioner is built
411: */
412: PetscErrorCode RHSMatrixLaplaciangllDM(TS ts, PetscReal t, Vec X, Mat A, Mat BB, PetscCtx ctx)
413: {
414: PetscReal **temp;
415: PetscReal vv;
416: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
417: PetscInt i, xs, xn, l, j;
418: PetscInt *rowsDM;
419: PetscBool flg = PETSC_FALSE;
421: PetscFunctionBeginUser;
422: PetscCall(PetscOptionsGetBool(NULL, NULL, "-gll_mf", &flg, NULL));
424: if (!flg) {
425: /*
426: Creates the element stiffness matrix for the given gll
427: */
428: PetscCall(PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
429: /* workaround for clang analyzer warning: Division by zero */
430: PetscCheck(appctx->param.N > 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Spectral element order should be > 1");
432: /* scale by the size of the element */
433: for (i = 0; i < appctx->param.N; i++) {
434: vv = -appctx->param.mu * 2.0 / appctx->param.Le;
435: for (j = 0; j < appctx->param.N; j++) temp[i][j] = temp[i][j] * vv;
436: }
438: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
439: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
441: xs = xs / (appctx->param.N - 1);
442: xn = xn / (appctx->param.N - 1);
444: PetscCall(PetscMalloc1(appctx->param.N, &rowsDM));
445: /*
446: loop over local elements
447: */
448: for (j = xs; j < xs + xn; j++) {
449: for (l = 0; l < appctx->param.N; l++) rowsDM[l] = 1 + (j - xs) * (appctx->param.N - 1) + l;
450: PetscCall(MatSetValuesLocal(A, appctx->param.N, rowsDM, appctx->param.N, rowsDM, &temp[0][0], ADD_VALUES));
451: }
452: PetscCall(PetscFree(rowsDM));
453: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
454: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
455: PetscCall(VecReciprocal(appctx->SEMop.mass));
456: PetscCall(MatDiagonalScale(A, appctx->SEMop.mass, NULL));
457: PetscCall(VecReciprocal(appctx->SEMop.mass));
459: PetscCall(PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
460: } else {
461: PetscCall(MatSetType(A, MATSHELL));
462: PetscCall(MatSetUp(A));
463: PetscCall(MatShellSetContext(A, appctx));
464: PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Laplacian));
465: }
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*
470: RHSMatrixAdvection - User-provided routine to compute the right-hand-side
471: matrix for the Advection (gradient) operator.
473: Input Parameters:
474: ts - the TS context
475: t - current time
476: global_in - global input vector
477: dummy - optional user-defined context, as set by TSetRHSJacobian()
479: Output Parameters:
480: AA - Jacobian matrix
481: BB - optionally different matrix used to construct the preconditioner
483: */
484: PetscErrorCode RHSMatrixAdvectiongllDM(TS ts, PetscReal t, Vec X, Mat A, Mat BB, PetscCtx ctx)
485: {
486: PetscReal **temp;
487: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
488: PetscInt xs, xn, l, j;
489: PetscInt *rowsDM;
490: PetscBool flg = PETSC_FALSE;
492: PetscFunctionBeginUser;
493: PetscCall(PetscOptionsGetBool(NULL, NULL, "-gll_mf", &flg, NULL));
495: if (!flg) {
496: /*
497: Creates the advection matrix for the given gll
498: */
499: PetscCall(PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
500: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
501: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
502: xs = xs / (appctx->param.N - 1);
503: xn = xn / (appctx->param.N - 1);
505: PetscCall(PetscMalloc1(appctx->param.N, &rowsDM));
506: for (j = xs; j < xs + xn; j++) {
507: for (l = 0; l < appctx->param.N; l++) rowsDM[l] = 1 + (j - xs) * (appctx->param.N - 1) + l;
508: PetscCall(MatSetValuesLocal(A, appctx->param.N, rowsDM, appctx->param.N, rowsDM, &temp[0][0], ADD_VALUES));
509: }
510: PetscCall(PetscFree(rowsDM));
511: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
512: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
514: PetscCall(VecReciprocal(appctx->SEMop.mass));
515: PetscCall(MatDiagonalScale(A, appctx->SEMop.mass, NULL));
516: PetscCall(VecReciprocal(appctx->SEMop.mass));
518: PetscCall(PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
519: } else {
520: PetscCall(MatSetType(A, MATSHELL));
521: PetscCall(MatSetUp(A));
522: PetscCall(MatShellSetContext(A, appctx));
523: PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Advection));
524: }
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: /*TEST
530: build:
531: requires: !complex
533: test:
534: suffix: 1
535: requires: !single
536: output_file: output/empty.out
538: test:
539: suffix: 2
540: nsize: 5
541: requires: !single
542: output_file: output/empty.out
544: test:
545: suffix: 3
546: requires: !single
547: args: -ts_view -ts_type beuler -gll_mf -pc_type none -ts_max_steps 5 -ts_monitor_error
549: test:
550: suffix: 4
551: requires: !single
552: args: -ts_view -ts_type beuler -pc_type none -ts_max_steps 5 -ts_monitor_error
554: TEST*/