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 */
100: PetscInitialize(&argc,&argv,(char*)0,help);
102: /* Create 1D DMStag for the solution, and set up. Note that you can supply many
103: command line options (see the man page for DMStagCreate1d)
104: */
105: DMStagCreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,1,1,DMSTAG_STENCIL_BOX,1,NULL,&dmSol);
106: DMSetFromOptions(dmSol);
107: DMSetUp(dmSol);
109: /* Create uniform coordinates. Note that in higher-dimensional examples,
110: coordinates are created differently.*/
111: DMStagSetUniformCoordinatesExplicit(dmSol,0.0,domainSize,0.0,0.0,0.0,0.0);
113: /* Determine boundary type */
114: DMStagGetBoundaryTypes(dmSol,&boundary,NULL,NULL);
116: /* Process command line options (depends on DMStag setup) */
117: a = 1.0; b = 2.0; c = 1.0;
118: PetscOptionsGetScalar(NULL,NULL,"-a",&a,NULL);
119: if (boundary == DM_BOUNDARY_PERIODIC) {
120: b = a;
121: c = 0.0;
122: } else {
123: PetscOptionsGetScalar(NULL,NULL,"-b",&b,NULL);
124: PetscOptionsGetScalar(NULL,NULL,"-c",&c,NULL);
125: }
127: /* Compute reference solution on the grid, using direct array access */
128: DMCreateGlobalVector(dmSol,&solRef);
129: DMGetLocalVector(dmSol,&solRefLocal);
130: DMStagVecGetArray(dmSol,solRefLocal,&arrSol);
131: DMGetCoordinateDM(dmSol,&dmCoordSol);
132: DMGetCoordinatesLocal(dmSol,&coordSolLocal);
133: DMStagVecGetArrayRead(dmCoordSol,coordSolLocal,&arrCoordSol);
134: DMStagGetCorners(dmSol,&start,NULL,NULL,&n,NULL,NULL,&nExtra,NULL,NULL);
136: /* Get the correct entries for each of our variables in local element-wise storage */
137: DMStagGetLocationSlot(dmSol,LEFT, 0,&iu);
138: DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip);
139: DMStagGetLocationSlot(dmCoordSol,LEFT, 0,&ixu);
140: DMStagGetLocationSlot(dmCoordSol,ELEMENT,0,&ixp);
141: for (e=start; e<start + n + nExtra; ++e) {
142: {
143: const PetscScalar coordu = arrCoordSol[e][ixu];
144: arrSol[e][iu] = a + (b - a - (c/2.0)) * coordu + (c/2.0)*coordu*coordu;
145: }
146: if (e < start+n) {
147: const PetscScalar coordp = arrCoordSol[e][ixp];
148: arrSol[e][ip] = b - a - (c/2.0) + c * coordp;
149: }
150: }
151: DMStagVecRestoreArrayRead(dmCoordSol,coordSolLocal,&arrCoordSol);
152: DMStagVecRestoreArray(dmSol,solRefLocal,&arrSol);
153: DMLocalToGlobal(dmSol,solRefLocal,INSERT_VALUES,solRef);
154: DMRestoreLocalVector(dmSol,&solRefLocal);
156: /* Create another 1D DMStag for the forcing term, and populate a field on it.
157: Here this is not really necessary, but in other contexts we may have auxiliary
158: fields which we use to construct the linear system.
160: This second DM represents the same physical domain, but has a different
161: "default section" (though the current implementation does NOT actually use
162: PetscSection). Since it is created as a derivative of the original DMStag,
163: we can be confident that it is compatible. One could check with DMGetCompatiblity() */
164: DMStagCreateCompatibleDMStag(dmSol,1,0,0,0,&dmForcing);
165: DMCreateGlobalVector(dmForcing,&f);
166: VecSet(f,c); /* Dummy for logic which depends on auxiliary data */
168: /* Assemble System */
169: DMCreateMatrix(dmSol,&A);
170: DMCreateGlobalVector(dmSol,&rhs);
171: DMGetLocalVector(dmForcing,&fLocal);
172: DMGlobalToLocal(dmForcing,f,INSERT_VALUES,fLocal);
174: /* Note: if iterating over all the elements, you will usually need to do something
175: special at one of the boundaries. You can either make use of the existence
176: of a "extra" partial dummy element on the right/top/front, or you can use a different stencil.
177: The construction of the reference solution above uses the first method,
178: so here we will use the second */
180: DMStagGetIsLastRank(dmSol,&isLastRank,NULL,NULL);
181: DMStagGetIsFirstRank(dmSol,&isFirstRank,NULL,NULL);
182: for (e = start; e<start+n; ++e) {
183: DMStagStencil pos[3];
184: PetscScalar val[3];
185: PetscInt idxLoc;
187: idxLoc = 0;
188: pos[idxLoc].i = e; /* This element in the 1d ordering */
189: pos[idxLoc].loc = ELEMENT; /* Element-centered dofs (p) */
190: pos[idxLoc].c = 0; /* Component 0 : first (and only) p dof */
191: val[idxLoc] = 0.0; /* p - u'(x) = 0 */
192: ++idxLoc;
194: if (isFirstRank && e == start) {
195: /* Special case on left boundary */
196: pos[idxLoc].i = e; /* This element in the 1d ordering */
197: pos[idxLoc].loc = LEFT; /* Left vertex */
198: pos[idxLoc].c = 0;
199: val[idxLoc] = a; /* u(0) = a */
200: ++idxLoc;
201: } else {
202: PetscScalar fVal;
203: /* Usual case - deal with velocity on left side of cell
204: Here, we obtain a value of f (even though it's constant here,
205: this demonstrates the more-realistic case of a pre-computed coefficient) */
206: pos[idxLoc].i = e; /* This element in the 1d ordering */
207: pos[idxLoc].loc = LEFT; /* vertex-centered dof (u) */
208: pos[idxLoc].c = 0;
210: DMStagVecGetValuesStencil(dmForcing,fLocal,1,&pos[idxLoc],&fVal);
212: val[idxLoc] = fVal; /* p'(x) = f, in interior */
213: ++idxLoc;
214: }
215: if (boundary != DM_BOUNDARY_PERIODIC && isLastRank && e == start+n-1) {
216: /* Special case on right boundary (in addition to usual case) */
217: pos[idxLoc].i = e; /* This element in the 1d ordering */
218: pos[idxLoc].loc = RIGHT;
219: pos[idxLoc].c = 0;
220: val[idxLoc] = b; /* u(1) = b */
221: ++idxLoc;
222: }
223: DMStagVecSetValuesStencil(dmSol,rhs,idxLoc,pos,val,INSERT_VALUES);
224: }
225: DMRestoreLocalVector(dmForcing,&fLocal);
226: VecAssemblyBegin(rhs);
227: VecAssemblyEnd(rhs);
229: /* Note: normally it would be more efficient to assemble the RHS and the matrix
230: in the same loop over elements, but we separate them for clarity here */
231: DMGetCoordinatesLocal(dmSol,&coordSolLocal);
232: for (e = start; e<start+n; ++e) {
234: /* Velocity is either a BC or an interior point */
235: if (isFirstRank && e == start) {
236: DMStagStencil row;
237: PetscScalar val;
239: row.i = e;
240: row.loc = LEFT;
241: row.c = 0;
242: val = 1.0;
243: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&val,INSERT_VALUES);
244: } else {
245: DMStagStencil row,col[3];
246: PetscScalar val[3],xp[2];
248: row.i = e;
249: row.loc = LEFT; /* In general, opt for LEFT/DOWN/BACK and iterate over elements */
250: row.c = 0;
252: col[0].i = e;
253: col[0].loc = ELEMENT;
254: col[0].c = 0;
256: col[1].i = e-1;
257: col[1].loc = ELEMENT;
258: col[1].c = 0;
260: DMStagVecGetValuesStencil(dmCoordSol,coordSolLocal,2,col,xp);
261: h = xp[0]- xp[1];
262: if (boundary == DM_BOUNDARY_PERIODIC && PetscRealPart(h) < 0.0) h += domainSize;
264: val[0] = 1.0/h;
265: val[1] = -1.0/h;
267: /* For convenience, we add an explicit 0 on the diagonal. This is a waste,
268: but it allows for easier use of a direct solver, if desired */
269: col[2].i = e;
270: col[2].loc = LEFT;
271: col[2].c = 0;
272: val[2] = 0.0;
274: DMStagMatSetValuesStencil(dmSol,A,1,&row,3,col,val,INSERT_VALUES);
275: }
277: /* Additional velocity point (BC) on the right */
278: if (boundary != DM_BOUNDARY_PERIODIC && isLastRank && e == start+n-1) {
279: DMStagStencil row;
280: PetscScalar val;
282: row.i = e;
283: row.loc = RIGHT;
284: row.c = 0;
285: val = 1.0;
286: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&val,INSERT_VALUES);
287: }
289: /* Equation on pressure (element) variables */
290: {
291: DMStagStencil row,col[3];
292: PetscScalar val[3],xu[2];
294: row.i = e;
295: row.loc = ELEMENT;
296: row.c = 0;
298: col[0].i = e;
299: col[0].loc = RIGHT;
300: col[0].c = 0;
302: col[1].i = e;
303: col[1].loc = LEFT;
304: col[1].c = 0;
306: DMStagVecGetValuesStencil(dmCoordSol,coordSolLocal,2,col,xu);
307: h = xu[0]- xu[1];
308: if (boundary == DM_BOUNDARY_PERIODIC && PetscRealPart(h) < 0.0) h += domainSize;
310: val[0] = -1.0/h;
311: val[1] = 1.0/h;
313: col[2].i = e;
314: col[2].loc = ELEMENT;
315: col[2].c = 0;
316: val[2] = 1.0;
318: DMStagMatSetValuesStencil(dmSol,A,1,&row,3,col,val,INSERT_VALUES);
319: }
320: }
321: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
322: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
324: /* Solve */
325: DMCreateGlobalVector(dmSol,&sol);
326: KSPCreate(PETSC_COMM_WORLD,&ksp);
327: KSPSetOperators(ksp,A,A);
328: KSPGetPC(ksp,&pc);
329: PCSetType(pc,PCJACOBI); /* A simple, but non-scalable, solver choice */
330: KSPSetFromOptions(ksp);
331: KSPSolve(ksp,rhs,sol);
333: /* View the components of the solution, demonstrating DMStagMigrateVec() */
334: {
335: DM dmVertsOnly,dmElementsOnly;
336: Vec u,p;
338: DMStagCreateCompatibleDMStag(dmSol,1,0,0,0,&dmVertsOnly);
339: DMStagCreateCompatibleDMStag(dmSol,0,1,0,0,&dmElementsOnly);
340: DMGetGlobalVector(dmVertsOnly,&u);
341: DMGetGlobalVector(dmElementsOnly,&p);
343: DMStagMigrateVec(dmSol,sol,dmVertsOnly,u);
344: DMStagMigrateVec(dmSol,sol,dmElementsOnly,p);
346: PetscObjectSetName((PetscObject)u,"Sol_u");
347: VecView(u,PETSC_VIEWER_STDOUT_WORLD);
348: PetscObjectSetName((PetscObject)p,"Sol_p");
349: VecView(p,PETSC_VIEWER_STDOUT_WORLD);
351: DMRestoreGlobalVector(dmVertsOnly,&u);
352: DMRestoreGlobalVector(dmElementsOnly,&p);
353: DMDestroy(&dmVertsOnly);
354: DMDestroy(&dmElementsOnly);
355: }
357: /* Check Solution */
358: {
359: Vec diff;
360: PetscReal normsolRef,errAbs,errRel;
362: VecDuplicate(sol,&diff);
363: VecCopy(sol,diff);
364: VecAXPY(diff,-1.0,solRef);
365: VecNorm(diff,NORM_2,&errAbs);
366: VecNorm(solRef,NORM_2,&normsolRef);
367: errRel = errAbs/normsolRef;
368: PetscPrintf(PETSC_COMM_WORLD,"Error (abs): %g\nError (rel): %g\n",(double)errAbs,(double)errRel);
369: VecDestroy(&diff);
370: }
372: /* Clean up and finalize PETSc */
373: KSPDestroy(&ksp);
374: VecDestroy(&sol);
375: VecDestroy(&solRef);
376: VecDestroy(&rhs);
377: VecDestroy(&f);
378: MatDestroy(&A);
379: DMDestroy(&dmSol);
380: DMDestroy(&dmForcing);
381: PetscFinalize();
382: return 0;
383: }
385: /*TEST
387: test:
388: suffix: 1
389: nsize: 7
390: args: -dm_view -stag_grid_x 11 -stag_stencil_type star -a 1.33 -b 7.22 -c 347.2 -ksp_monitor_short
392: test:
393: suffix: periodic
394: nsize: 3
395: args: -dm_view -stag_grid_x 13 -stag_boundary_type_x periodic -a 1.1234
397: test:
398: suffix: periodic_seq
399: nsize: 1
400: args: -dm_view -stag_grid_x 13 -stag_boundary_type_x periodic -a 1.1234
402: test:
403: suffix: ghosted_vacuous
404: nsize: 3
405: args: -dm_view -stag_grid_x 13 -stag_boundary_type_x ghosted -a 1.1234
407: TEST*/