Actual source code: ex1.c
1: static char help[] = "Solve a toy 1D problem on a staggered grid.\n\
2: Accepts command line options -a, -b, and -c\n\
3: and approximately solves\n\
4: u''(x) = c, u(0) = 1, u(1) = b\n\n";
5: /*
7: To demonstrate the basic functionality of DMStag, solves a second-order ODE,
9: u''(x) = f(x), 0 < x < 1
10: u(0) = a
11: u(1) = b
13: in mixed form, that is by introducing an auxiliary variable p
15: p'(x) = f(x), p - u'(x) = 0, 0 < x < 1
16: u(0) = a
17: u(1) = b
19: For f == c, the solution is
21: u(x) = a + (b - a - (c/2)) x + (c/2) x^2
22: p(x) = (b - a - (c/2)) + c x
24: To find an approximate solution, discretize by storing values of p in
25: elements and values of u on their boundaries, and using first-order finite
26: differences.
28: This should in fact produce a (nodal) solution with no discretization error,
29: so differences from the reference solution will be restricted to those induced
30: by floating point operations (in particular, the finite differences) and the
31: accuracy of the linear solve.
33: Parameters for the main grid can be controlled with command line options, e.g.
35: -stag_grid_x 10
37: In particular to notice in this example are the two methods of indexing. The
38: first is analogous to the use of MatStencil with DMDA, and the second is
39: analogous to the use of DMDAVecGetArrayDOF().
41: The first, recommended for ease of use, is based on naming an element in the
42: global grid, a location in its support, and a component. For example,
43: one might consider element e, the left side (a vertex in 1d), and the first
44: component (index 0). This is accomplished by populating a DMStagStencil struct,
45: e.g.
47: DMStagStencil stencil;
48: stencil.i = i
49: stencil.loc = DMSTAG_LEFT;
50: stencil.c = 0
52: Note that below, for convenenience, we #define an alias LEFT for DMSTAG_LEFT.
54: The second, which ensures maximum efficiency, makes use of the underlying
55: block structure of local DMStag-derived vectors, and requires the user to
56: obtain the correct local offset for the degrees of freedom they would like to
57: use. This is made simple with the helper function DMStagGetLocationSlot().
59: Note that the linear system being solved is indefinite, so is not entirely
60: trivial to invert. The default solver here (GMRES/Jacobi) is a poor choice,
61: made to avoid depending on an external package. To solve a larger system,
62: the usual method for a 1-d problem such as this is to choose a sophisticated
63: direct solver, e.g. configure --download-suitesparse and run
65: $PETSC_DIR/$PETSC_ARCH/bin/mpiexec -n 3 ./stag_ex2 -stag_grid_x 100 -pc_type lu -pc_factor_mat_solver_package umfpack
67: You can also impose a periodic boundary condition, in which case -b and -c are
68: ignored; b = a and c = 0.0 are used, giving a constant u == a , p == 0.
70: -stag_boundary_type_x periodic
72: */
73: #include <petscdm.h>
74: #include <petscksp.h>
75: #include <petscdmstag.h>
77: /* Shorter, more convenient names for DMStagStencilLocation entries */
78: #define LEFT DMSTAG_LEFT
79: #define RIGHT DMSTAG_RIGHT
80: #define ELEMENT DMSTAG_ELEMENT
82: int main(int argc, char **argv)
83: {
84: DM dmSol, dmForcing;
85: DM dmCoordSol;
86: Vec sol, solRef, solRefLocal, f, fLocal, rhs, coordSolLocal;
87: Mat A;
88: PetscScalar a, b, c, h;
89: KSP ksp;
90: PC pc;
91: PetscInt start, n, e, nExtra;
92: PetscInt iu, ip, ixu, ixp;
93: PetscBool isLastRank, isFirstRank;
94: PetscScalar **arrSol, **arrCoordSol;
95: DMBoundaryType boundary;
97: const PetscReal domainSize = 1.0;
99: /* Initialize PETSc */
101: PetscInitialize(&argc, &argv, (char *)0, help);
103: /* Create 1D DMStag for the solution, and set up. Note that you can supply many
104: command line options (see the man page for DMStagCreate1d)
105: */
106: DMStagCreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 3, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, &dmSol);
107: DMSetFromOptions(dmSol);
108: DMSetUp(dmSol);
110: /* Create uniform coordinates. Note that in higher-dimensional examples,
111: coordinates are created differently.*/
112: DMStagSetUniformCoordinatesExplicit(dmSol, 0.0, domainSize, 0.0, 0.0, 0.0, 0.0);
114: /* Determine boundary type */
115: DMStagGetBoundaryTypes(dmSol, &boundary, NULL, NULL);
117: /* Process command line options (depends on DMStag setup) */
118: a = 1.0;
119: b = 2.0;
120: c = 1.0;
121: PetscOptionsGetScalar(NULL, NULL, "-a", &a, NULL);
122: if (boundary == DM_BOUNDARY_PERIODIC) {
123: b = a;
124: c = 0.0;
125: } else {
126: PetscOptionsGetScalar(NULL, NULL, "-b", &b, NULL);
127: PetscOptionsGetScalar(NULL, NULL, "-c", &c, NULL);
128: }
130: /* Compute reference solution on the grid, using direct array access */
131: DMCreateGlobalVector(dmSol, &solRef);
132: DMGetLocalVector(dmSol, &solRefLocal);
133: DMStagVecGetArray(dmSol, solRefLocal, &arrSol);
134: DMGetCoordinateDM(dmSol, &dmCoordSol);
135: DMGetCoordinatesLocal(dmSol, &coordSolLocal);
136: DMStagVecGetArrayRead(dmCoordSol, coordSolLocal, &arrCoordSol);
137: DMStagGetCorners(dmSol, &start, NULL, NULL, &n, NULL, NULL, &nExtra, NULL, NULL);
139: /* Get the correct entries for each of our variables in local element-wise storage */
140: DMStagGetLocationSlot(dmSol, LEFT, 0, &iu);
141: DMStagGetLocationSlot(dmSol, ELEMENT, 0, &ip);
142: DMStagGetLocationSlot(dmCoordSol, LEFT, 0, &ixu);
143: DMStagGetLocationSlot(dmCoordSol, ELEMENT, 0, &ixp);
144: for (e = start; e < start + n + nExtra; ++e) {
145: {
146: const PetscScalar coordu = arrCoordSol[e][ixu];
147: arrSol[e][iu] = a + (b - a - (c / 2.0)) * coordu + (c / 2.0) * coordu * coordu;
148: }
149: if (e < start + n) {
150: const PetscScalar coordp = arrCoordSol[e][ixp];
151: arrSol[e][ip] = b - a - (c / 2.0) + c * coordp;
152: }
153: }
154: DMStagVecRestoreArrayRead(dmCoordSol, coordSolLocal, &arrCoordSol);
155: DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol);
156: DMLocalToGlobal(dmSol, solRefLocal, INSERT_VALUES, solRef);
157: DMRestoreLocalVector(dmSol, &solRefLocal);
159: /* Create another 1D DMStag for the forcing term, and populate a field on it.
160: Here this is not really necessary, but in other contexts we may have auxiliary
161: fields which we use to construct the linear system.
163: This second DM represents the same physical domain, but has a different
164: "default section" (though the current implementation does NOT actually use
165: PetscSection). Since it is created as a derivative of the original DMStag,
166: we can be confident that it is compatible. One could check with DMGetCompatiblity() */
167: DMStagCreateCompatibleDMStag(dmSol, 1, 0, 0, 0, &dmForcing);
168: DMCreateGlobalVector(dmForcing, &f);
169: VecSet(f, c); /* Dummy for logic which depends on auxiliary data */
171: /* Assemble System */
172: DMCreateMatrix(dmSol, &A);
173: DMCreateGlobalVector(dmSol, &rhs);
174: DMGetLocalVector(dmForcing, &fLocal);
175: DMGlobalToLocal(dmForcing, f, INSERT_VALUES, fLocal);
177: /* Note: if iterating over all the elements, you will usually need to do something
178: special at one of the boundaries. You can either make use of the existence
179: of a "extra" partial dummy element on the right/top/front, or you can use a different stencil.
180: The construction of the reference solution above uses the first method,
181: so here we will use the second */
183: DMStagGetIsLastRank(dmSol, &isLastRank, NULL, NULL);
184: DMStagGetIsFirstRank(dmSol, &isFirstRank, NULL, NULL);
185: for (e = start; e < start + n; ++e) {
186: DMStagStencil pos[3];
187: PetscScalar val[3];
188: PetscInt idxLoc;
190: idxLoc = 0;
191: pos[idxLoc].i = e; /* This element in the 1d ordering */
192: pos[idxLoc].loc = ELEMENT; /* Element-centered dofs (p) */
193: pos[idxLoc].c = 0; /* Component 0 : first (and only) p dof */
194: val[idxLoc] = 0.0; /* p - u'(x) = 0 */
195: ++idxLoc;
197: if (isFirstRank && e == start) {
198: /* Special case on left boundary */
199: pos[idxLoc].i = e; /* This element in the 1d ordering */
200: pos[idxLoc].loc = LEFT; /* Left vertex */
201: pos[idxLoc].c = 0;
202: val[idxLoc] = a; /* u(0) = a */
203: ++idxLoc;
204: } else {
205: PetscScalar fVal;
206: /* Usual case - deal with velocity on left side of cell
207: Here, we obtain a value of f (even though it's constant here,
208: this demonstrates the more-realistic case of a pre-computed coefficient) */
209: pos[idxLoc].i = e; /* This element in the 1d ordering */
210: pos[idxLoc].loc = LEFT; /* vertex-centered dof (u) */
211: pos[idxLoc].c = 0;
213: DMStagVecGetValuesStencil(dmForcing, fLocal, 1, &pos[idxLoc], &fVal);
215: val[idxLoc] = fVal; /* p'(x) = f, in interior */
216: ++idxLoc;
217: }
218: if (boundary != DM_BOUNDARY_PERIODIC && isLastRank && e == start + n - 1) {
219: /* Special case on right boundary (in addition to usual case) */
220: pos[idxLoc].i = e; /* This element in the 1d ordering */
221: pos[idxLoc].loc = RIGHT;
222: pos[idxLoc].c = 0;
223: val[idxLoc] = b; /* u(1) = b */
224: ++idxLoc;
225: }
226: DMStagVecSetValuesStencil(dmSol, rhs, idxLoc, pos, val, INSERT_VALUES);
227: }
228: DMRestoreLocalVector(dmForcing, &fLocal);
229: VecAssemblyBegin(rhs);
230: VecAssemblyEnd(rhs);
232: /* Note: normally it would be more efficient to assemble the RHS and the matrix
233: in the same loop over elements, but we separate them for clarity here */
234: DMGetCoordinatesLocal(dmSol, &coordSolLocal);
235: for (e = start; e < start + n; ++e) {
236: /* Velocity is either a BC or an interior point */
237: if (isFirstRank && e == start) {
238: DMStagStencil row;
239: PetscScalar val;
241: row.i = e;
242: row.loc = LEFT;
243: row.c = 0;
244: val = 1.0;
245: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &val, INSERT_VALUES);
246: } else {
247: DMStagStencil row, col[3];
248: PetscScalar val[3], xp[2];
250: row.i = e;
251: row.loc = LEFT; /* In general, opt for LEFT/DOWN/BACK and iterate over elements */
252: row.c = 0;
254: col[0].i = e;
255: col[0].loc = ELEMENT;
256: col[0].c = 0;
258: col[1].i = e - 1;
259: col[1].loc = ELEMENT;
260: col[1].c = 0;
262: DMStagVecGetValuesStencil(dmCoordSol, coordSolLocal, 2, col, xp);
263: h = xp[0] - xp[1];
264: if (boundary == DM_BOUNDARY_PERIODIC && PetscRealPart(h) < 0.0) h += domainSize;
266: val[0] = 1.0 / h;
267: val[1] = -1.0 / h;
269: /* For convenience, we add an explicit 0 on the diagonal. This is a waste,
270: but it allows for easier use of a direct solver, if desired */
271: col[2].i = e;
272: col[2].loc = LEFT;
273: col[2].c = 0;
274: val[2] = 0.0;
276: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 3, col, val, INSERT_VALUES);
277: }
279: /* Additional velocity point (BC) on the right */
280: if (boundary != DM_BOUNDARY_PERIODIC && isLastRank && e == start + n - 1) {
281: DMStagStencil row;
282: PetscScalar val;
284: row.i = e;
285: row.loc = RIGHT;
286: row.c = 0;
287: val = 1.0;
288: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &val, INSERT_VALUES);
289: }
291: /* Equation on pressure (element) variables */
292: {
293: DMStagStencil row, col[3];
294: PetscScalar val[3], xu[2];
296: row.i = e;
297: row.loc = ELEMENT;
298: row.c = 0;
300: col[0].i = e;
301: col[0].loc = RIGHT;
302: col[0].c = 0;
304: col[1].i = e;
305: col[1].loc = LEFT;
306: col[1].c = 0;
308: DMStagVecGetValuesStencil(dmCoordSol, coordSolLocal, 2, col, xu);
309: h = xu[0] - xu[1];
310: if (boundary == DM_BOUNDARY_PERIODIC && PetscRealPart(h) < 0.0) h += domainSize;
312: val[0] = -1.0 / h;
313: val[1] = 1.0 / h;
315: col[2].i = e;
316: col[2].loc = ELEMENT;
317: col[2].c = 0;
318: val[2] = 1.0;
320: DMStagMatSetValuesStencil(dmSol, A, 1, &row, 3, col, val, INSERT_VALUES);
321: }
322: }
323: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
324: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
326: /* Solve */
327: DMCreateGlobalVector(dmSol, &sol);
328: KSPCreate(PETSC_COMM_WORLD, &ksp);
329: KSPSetOperators(ksp, A, A);
330: KSPGetPC(ksp, &pc);
331: PCSetType(pc, PCJACOBI); /* A simple, but non-scalable, solver choice */
332: KSPSetFromOptions(ksp);
333: KSPSolve(ksp, rhs, sol);
335: /* View the components of the solution, demonstrating DMStagMigrateVec() */
336: {
337: DM dmVertsOnly, dmElementsOnly;
338: Vec u, p;
340: DMStagCreateCompatibleDMStag(dmSol, 1, 0, 0, 0, &dmVertsOnly);
341: DMStagCreateCompatibleDMStag(dmSol, 0, 1, 0, 0, &dmElementsOnly);
342: DMGetGlobalVector(dmVertsOnly, &u);
343: DMGetGlobalVector(dmElementsOnly, &p);
345: DMStagMigrateVec(dmSol, sol, dmVertsOnly, u);
346: DMStagMigrateVec(dmSol, sol, dmElementsOnly, p);
348: PetscObjectSetName((PetscObject)u, "Sol_u");
349: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
350: PetscObjectSetName((PetscObject)p, "Sol_p");
351: VecView(p, PETSC_VIEWER_STDOUT_WORLD);
353: DMRestoreGlobalVector(dmVertsOnly, &u);
354: DMRestoreGlobalVector(dmElementsOnly, &p);
355: DMDestroy(&dmVertsOnly);
356: DMDestroy(&dmElementsOnly);
357: }
359: /* Check Solution */
360: {
361: Vec diff;
362: PetscReal normsolRef, errAbs, errRel;
364: VecDuplicate(sol, &diff);
365: VecCopy(sol, diff);
366: VecAXPY(diff, -1.0, solRef);
367: VecNorm(diff, NORM_2, &errAbs);
368: VecNorm(solRef, NORM_2, &normsolRef);
369: errRel = errAbs / normsolRef;
370: PetscPrintf(PETSC_COMM_WORLD, "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel);
371: VecDestroy(&diff);
372: }
374: /* Clean up and finalize PETSc */
375: KSPDestroy(&ksp);
376: VecDestroy(&sol);
377: VecDestroy(&solRef);
378: VecDestroy(&rhs);
379: VecDestroy(&f);
380: MatDestroy(&A);
381: DMDestroy(&dmSol);
382: DMDestroy(&dmForcing);
383: PetscFinalize();
384: return 0;
385: }
387: /*TEST
389: test:
390: suffix: 1
391: nsize: 7
392: args: -dm_view -stag_grid_x 11 -stag_stencil_type star -a 1.33 -b 7.22 -c 347.2 -ksp_monitor_short
394: test:
395: suffix: periodic
396: nsize: 3
397: args: -dm_view -stag_grid_x 13 -stag_boundary_type_x periodic -a 1.1234
399: test:
400: suffix: periodic_seq
401: nsize: 1
402: args: -dm_view -stag_grid_x 13 -stag_boundary_type_x periodic -a 1.1234
404: test:
405: suffix: ghosted_vacuous
406: nsize: 3
407: args: -dm_view -stag_grid_x 13 -stag_boundary_type_x ghosted -a 1.1234
409: TEST*/