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*/