Actual source code: ex13.c

  1: static char help[] = "Poisson Problem in 2d and 3d with finite elements.\n\
  2: We solve the Poisson problem in a rectangular\n\
  3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
  4: This example supports automatic convergence estimation\n\
  5: and eventually adaptivity.\n\n\n";

  7: #include <petscdmplex.h>
  8: #include <petscsnes.h>
  9: #include <petscds.h>
 10: #include <petscconvest.h>

 12: typedef struct {
 13:   /* Domain and mesh definition */
 14:   PetscBool spectral;    /* Look at the spectrum along planes in the solution */
 15:   PetscBool shear;       /* Shear the domain */
 16:   PetscBool adjoint;     /* Solve the adjoint problem */
 17:   PetscBool homogeneous; /* Use homogeneous boudnary conditions */
 18:   PetscBool viewError;   /* Output the solution error */
 19: } AppCtx;

 21: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 22: {
 23:   *u = 0.0;
 24:   return 0;
 25: }

 27: static PetscErrorCode trig_inhomogeneous_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 28: {
 29:   PetscInt d;
 30:   *u = 0.0;
 31:   for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]);
 32:   return 0;
 33: }

 35: static PetscErrorCode trig_homogeneous_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 36: {
 37:   PetscInt d;
 38:   *u = 1.0;
 39:   for (d = 0; d < dim; ++d) *u *= PetscSinReal(2.0 * PETSC_PI * x[d]);
 40:   return 0;
 41: }

 43: /* Compute integral of (residual of solution)*(adjoint solution - projection of adjoint solution) */
 44: static void obj_error_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
 45: {
 46:   obj[0] = a[aOff[0]] * (u[0] - a[aOff[1]]);
 47: }

 49: static void f0_trig_inhomogeneous_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 50: {
 51:   PetscInt d;
 52:   for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
 53: }

 55: static void f0_trig_homogeneous_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 56: {
 57:   PetscInt d;
 58:   for (d = 0; d < dim; ++d) {
 59:     PetscScalar v = 1.;
 60:     for (PetscInt e = 0; e < dim; e++) {
 61:       if (e == d) {
 62:         v *= -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
 63:       } else {
 64:         v *= PetscSinReal(2.0 * PETSC_PI * x[d]);
 65:       }
 66:     }
 67:     f0[0] += v;
 68:   }
 69: }

 71: static void f0_unity_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 72: {
 73:   f0[0] = 1.0;
 74: }

 76: static void f0_identityaux_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 77: {
 78:   f0[0] = a[0];
 79: }

 81: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
 82: {
 83:   PetscInt d;
 84:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
 85: }

 87: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
 88: {
 89:   PetscInt d;
 90:   for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
 91: }

 93: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 94: {
 96:   options->shear       = PETSC_FALSE;
 97:   options->spectral    = PETSC_FALSE;
 98:   options->adjoint     = PETSC_FALSE;
 99:   options->homogeneous = PETSC_FALSE;
100:   options->viewError   = PETSC_FALSE;

102:   PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
103:   PetscOptionsBool("-shear", "Shear the domain", "ex13.c", options->shear, &options->shear, NULL);
104:   PetscOptionsBool("-spectral", "Look at the spectrum along planes of the solution", "ex13.c", options->spectral, &options->spectral, NULL);
105:   PetscOptionsBool("-adjoint", "Solve the adjoint problem", "ex13.c", options->adjoint, &options->adjoint, NULL);
106:   PetscOptionsBool("-homogeneous", "Use homogeneous boundary conditions", "ex13.c", options->homogeneous, &options->homogeneous, NULL);
107:   PetscOptionsBool("-error_view", "Output the solution error", "ex13.c", options->viewError, &options->viewError, NULL);
108:   PetscOptionsEnd();
109:   return 0;
110: }

112: static PetscErrorCode CreateSpectralPlanes(DM dm, PetscInt numPlanes, const PetscInt planeDir[], const PetscReal planeCoord[], AppCtx *user)
113: {
114:   PetscSection       coordSection;
115:   Vec                coordinates;
116:   const PetscScalar *coords;
117:   PetscInt           dim, p, vStart, vEnd, v;

120:   DMGetCoordinateDim(dm, &dim);
121:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
122:   DMGetCoordinatesLocal(dm, &coordinates);
123:   DMGetCoordinateSection(dm, &coordSection);
124:   VecGetArrayRead(coordinates, &coords);
125:   for (p = 0; p < numPlanes; ++p) {
126:     DMLabel label;
127:     char    name[PETSC_MAX_PATH_LEN];

129:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p);
130:     DMCreateLabel(dm, name);
131:     DMGetLabel(dm, name, &label);
132:     DMLabelAddStratum(label, 1);
133:     for (v = vStart; v < vEnd; ++v) {
134:       PetscInt off;

136:       PetscSectionGetOffset(coordSection, v, &off);
137:       if (PetscAbsReal(planeCoord[p] - PetscRealPart(coords[off + planeDir[p]])) < PETSC_SMALL) DMLabelSetValue(label, v, 1);
138:     }
139:   }
140:   VecRestoreArrayRead(coordinates, &coords);
141:   return 0;
142: }

144: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
145: {
147:   DMCreate(comm, dm);
148:   DMSetType(*dm, DMPLEX);
149:   DMSetFromOptions(*dm);
150:   if (user->shear) DMPlexShearGeometry(*dm, DM_X, NULL);
151:   DMSetApplicationContext(*dm, user);
152:   DMViewFromOptions(*dm, NULL, "-dm_view");
153:   if (user->spectral) {
154:     PetscInt  planeDir[2]   = {0, 1};
155:     PetscReal planeCoord[2] = {0., 1.};

157:     CreateSpectralPlanes(*dm, 2, planeDir, planeCoord, user);
158:   }
159:   return 0;
160: }

162: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
163: {
164:   PetscDS        ds;
165:   DMLabel        label;
166:   const PetscInt id                                                                             = 1;
167:   PetscPointFunc f0                                                                             = user->homogeneous ? f0_trig_homogeneous_u : f0_trig_inhomogeneous_u;
168:   PetscErrorCode (*ex)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = user->homogeneous ? trig_homogeneous_u : trig_inhomogeneous_u;

171:   DMGetDS(dm, &ds);
172:   PetscDSSetResidual(ds, 0, f0, f1_u);
173:   PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);
174:   PetscDSSetExactSolution(ds, 0, ex, user);
175:   DMGetLabel(dm, "marker", &label);
176:   if (label) DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))ex, NULL, user, NULL);
177:   return 0;
178: }

180: static PetscErrorCode SetupAdjointProblem(DM dm, AppCtx *user)
181: {
182:   PetscDS        ds;
183:   DMLabel        label;
184:   const PetscInt id = 1;

187:   DMGetDS(dm, &ds);
188:   PetscDSSetResidual(ds, 0, f0_unity_u, f1_u);
189:   PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);
190:   PetscDSSetObjective(ds, 0, obj_error_u);
191:   DMGetLabel(dm, "marker", &label);
192:   DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))zero, NULL, user, NULL);
193:   return 0;
194: }

196: static PetscErrorCode SetupErrorProblem(DM dm, AppCtx *user)
197: {
198:   PetscDS prob;

201:   DMGetDS(dm, &prob);
202:   return 0;
203: }

205: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
206: {
207:   DM             cdm = dm;
208:   PetscFE        fe;
209:   DMPolytopeType ct;
210:   PetscBool      simplex;
211:   PetscInt       dim, cStart;
212:   char           prefix[PETSC_MAX_PATH_LEN];

215:   DMGetDimension(dm, &dim);
216:   DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
217:   DMPlexGetCellType(dm, cStart, &ct);
218:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
219:   /* Create finite element */
220:   PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
221:   PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe);
222:   PetscObjectSetName((PetscObject)fe, name);
223:   /* Set discretization and boundary conditions for each mesh */
224:   DMSetField(dm, 0, NULL, (PetscObject)fe);
225:   DMCreateDS(dm);
226:   (*setup)(dm, user);
227:   while (cdm) {
228:     DMCopyDisc(dm, cdm);
229:     /* TODO: Check whether the boundary of coarse meshes is marked */
230:     DMGetCoarseDM(cdm, &cdm);
231:   }
232:   PetscFEDestroy(&fe);
233:   return 0;
234: }

236: static PetscErrorCode ComputeSpectral(DM dm, Vec u, PetscInt numPlanes, const PetscInt planeDir[], const PetscReal planeCoord[], AppCtx *user)
237: {
238:   MPI_Comm           comm;
239:   PetscSection       coordSection, section;
240:   Vec                coordinates, uloc;
241:   const PetscScalar *coords, *array;
242:   PetscInt           p;
243:   PetscMPIInt        size, rank;

246:   PetscObjectGetComm((PetscObject)dm, &comm);
247:   MPI_Comm_size(comm, &size);
248:   MPI_Comm_rank(comm, &rank);
249:   DMGetLocalVector(dm, &uloc);
250:   DMGlobalToLocalBegin(dm, u, INSERT_VALUES, uloc);
251:   DMGlobalToLocalEnd(dm, u, INSERT_VALUES, uloc);
252:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, uloc, 0.0, NULL, NULL, NULL);
253:   VecViewFromOptions(uloc, NULL, "-sol_view");
254:   DMGetLocalSection(dm, &section);
255:   VecGetArrayRead(uloc, &array);
256:   DMGetCoordinatesLocal(dm, &coordinates);
257:   DMGetCoordinateSection(dm, &coordSection);
258:   VecGetArrayRead(coordinates, &coords);
259:   for (p = 0; p < numPlanes; ++p) {
260:     DMLabel         label;
261:     char            name[PETSC_MAX_PATH_LEN];
262:     Mat             F;
263:     Vec             x, y;
264:     IS              stratum;
265:     PetscReal      *ray, *gray;
266:     PetscScalar    *rvals, *svals, *gsvals;
267:     PetscInt       *perm, *nperm;
268:     PetscInt        n, N, i, j, off, offu;
269:     const PetscInt *points;

271:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p);
272:     DMGetLabel(dm, name, &label);
273:     DMLabelGetStratumIS(label, 1, &stratum);
274:     ISGetLocalSize(stratum, &n);
275:     ISGetIndices(stratum, &points);
276:     PetscMalloc2(n, &ray, n, &svals);
277:     for (i = 0; i < n; ++i) {
278:       PetscSectionGetOffset(coordSection, points[i], &off);
279:       PetscSectionGetOffset(section, points[i], &offu);
280:       ray[i]   = PetscRealPart(coords[off + ((planeDir[p] + 1) % 2)]);
281:       svals[i] = array[offu];
282:     }
283:     /* Gather the ray data to proc 0 */
284:     if (size > 1) {
285:       PetscMPIInt *cnt, *displs, p;

287:       PetscCalloc2(size, &cnt, size, &displs);
288:       MPI_Gather(&n, 1, MPIU_INT, cnt, 1, MPIU_INT, 0, comm);
289:       for (p = 1; p < size; ++p) displs[p] = displs[p - 1] + cnt[p - 1];
290:       N = displs[size - 1] + cnt[size - 1];
291:       PetscMalloc2(N, &gray, N, &gsvals);
292:       MPI_Gatherv(ray, n, MPIU_REAL, gray, cnt, displs, MPIU_REAL, 0, comm);
293:       MPI_Gatherv(svals, n, MPIU_SCALAR, gsvals, cnt, displs, MPIU_SCALAR, 0, comm);
294:       PetscFree2(cnt, displs);
295:     } else {
296:       N      = n;
297:       gray   = ray;
298:       gsvals = svals;
299:     }
300:     if (rank == 0) {
301:       /* Sort point along ray */
302:       PetscMalloc2(N, &perm, N, &nperm);
303:       for (i = 0; i < N; ++i) perm[i] = i;
304:       PetscSortRealWithPermutation(N, gray, perm);
305:       /* Count duplicates and squish mapping */
306:       nperm[0] = perm[0];
307:       for (i = 1, j = 1; i < N; ++i) {
308:         if (PetscAbsReal(gray[perm[i]] - gray[perm[i - 1]]) > PETSC_SMALL) nperm[j++] = perm[i];
309:       }
310:       /* Create FFT structs */
311:       MatCreateFFT(PETSC_COMM_SELF, 1, &j, MATFFTW, &F);
312:       MatCreateVecs(F, &x, &y);
313:       PetscObjectSetName((PetscObject)y, name);
314:       VecGetArray(x, &rvals);
315:       for (i = 0, j = 0; i < N; ++i) {
316:         if (i > 0 && PetscAbsReal(gray[perm[i]] - gray[perm[i - 1]]) < PETSC_SMALL) continue;
317:         rvals[j] = gsvals[nperm[j]];
318:         ++j;
319:       }
320:       PetscFree2(perm, nperm);
321:       if (size > 1) PetscFree2(gray, gsvals);
322:       VecRestoreArray(x, &rvals);
323:       /* Do FFT along the ray */
324:       MatMult(F, x, y);
325:       /* Chop FFT */
326:       VecChop(y, PETSC_SMALL);
327:       VecViewFromOptions(x, NULL, "-real_view");
328:       VecViewFromOptions(y, NULL, "-fft_view");
329:       VecDestroy(&x);
330:       VecDestroy(&y);
331:       MatDestroy(&F);
332:     }
333:     ISRestoreIndices(stratum, &points);
334:     ISDestroy(&stratum);
335:     PetscFree2(ray, svals);
336:   }
337:   VecRestoreArrayRead(coordinates, &coords);
338:   VecRestoreArrayRead(uloc, &array);
339:   DMRestoreLocalVector(dm, &uloc);
340:   return 0;
341: }

343: int main(int argc, char **argv)
344: {
345:   DM     dm;   /* Problem specification */
346:   SNES   snes; /* Nonlinear solver */
347:   Vec    u;    /* Solutions */
348:   AppCtx user; /* User-defined work context */

351:   PetscInitialize(&argc, &argv, NULL, help);
352:   ProcessOptions(PETSC_COMM_WORLD, &user);
353:   /* Primal system */
354:   SNESCreate(PETSC_COMM_WORLD, &snes);
355:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
356:   SNESSetDM(snes, dm);
357:   SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);
358:   DMCreateGlobalVector(dm, &u);
359:   VecSet(u, 0.0);
360:   PetscObjectSetName((PetscObject)u, "potential");
361:   DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
362:   SNESSetFromOptions(snes);
363:   SNESSolve(snes, NULL, u);
364:   SNESGetSolution(snes, &u);
365:   VecViewFromOptions(u, NULL, "-potential_view");
366:   if (user.viewError) {
367:     PetscErrorCode (*sol)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
368:     void     *ctx;
369:     PetscDS   ds;
370:     PetscReal error;
371:     PetscInt  N;

373:     DMGetDS(dm, &ds);
374:     PetscDSGetExactSolution(ds, 0, &sol, &ctx);
375:     VecGetSize(u, &N);
376:     DMComputeL2Diff(dm, 0.0, &sol, &ctx, u, &error);
377:     PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " L2 error: %g\n", N, (double)error);
378:   }
379:   if (user.spectral) {
380:     PetscInt  planeDir[2]   = {0, 1};
381:     PetscReal planeCoord[2] = {0., 1.};

383:     ComputeSpectral(dm, u, 2, planeDir, planeCoord, &user);
384:   }
385:   /* Adjoint system */
386:   if (user.adjoint) {
387:     DM   dmAdj;
388:     SNES snesAdj;
389:     Vec  uAdj;

391:     SNESCreate(PETSC_COMM_WORLD, &snesAdj);
392:     PetscObjectSetOptionsPrefix((PetscObject)snesAdj, "adjoint_");
393:     DMClone(dm, &dmAdj);
394:     SNESSetDM(snesAdj, dmAdj);
395:     SetupDiscretization(dmAdj, "adjoint", SetupAdjointProblem, &user);
396:     DMCreateGlobalVector(dmAdj, &uAdj);
397:     VecSet(uAdj, 0.0);
398:     PetscObjectSetName((PetscObject)uAdj, "adjoint");
399:     DMPlexSetSNESLocalFEM(dmAdj, &user, &user, &user);
400:     SNESSetFromOptions(snesAdj);
401:     SNESSolve(snesAdj, NULL, uAdj);
402:     SNESGetSolution(snesAdj, &uAdj);
403:     VecViewFromOptions(uAdj, NULL, "-adjoint_view");
404:     /* Error representation */
405:     {
406:       DM        dmErr, dmErrAux, dms[2];
407:       Vec       errorEst, errorL2, uErr, uErrLoc, uAdjLoc, uAdjProj;
408:       IS       *subis;
409:       PetscReal errorEstTot, errorL2Norm, errorL2Tot;
410:       PetscInt  N, i;
411:       PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {user.homogeneous ? trig_homogeneous_u : trig_inhomogeneous_u};
412:       void (*identity[1])(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]) = {f0_identityaux_u};
413:       void *ctxs[1] = {0};

415:       ctxs[0] = &user;
416:       DMClone(dm, &dmErr);
417:       SetupDiscretization(dmErr, "error", SetupErrorProblem, &user);
418:       DMGetGlobalVector(dmErr, &errorEst);
419:       DMGetGlobalVector(dmErr, &errorL2);
420:       /*   Compute auxiliary data (solution and projection of adjoint solution) */
421:       DMGetLocalVector(dmAdj, &uAdjLoc);
422:       DMGlobalToLocalBegin(dmAdj, uAdj, INSERT_VALUES, uAdjLoc);
423:       DMGlobalToLocalEnd(dmAdj, uAdj, INSERT_VALUES, uAdjLoc);
424:       DMGetGlobalVector(dm, &uAdjProj);
425:       DMSetAuxiliaryVec(dm, NULL, 0, 0, uAdjLoc);
426:       DMProjectField(dm, 0.0, u, identity, INSERT_VALUES, uAdjProj);
427:       DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL);
428:       DMRestoreLocalVector(dmAdj, &uAdjLoc);
429:       /*   Attach auxiliary data */
430:       dms[0] = dm;
431:       dms[1] = dm;
432:       DMCreateSuperDM(dms, 2, &subis, &dmErrAux);
433:       if (0) {
434:         PetscSection sec;

436:         DMGetLocalSection(dms[0], &sec);
437:         PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD);
438:         DMGetLocalSection(dms[1], &sec);
439:         PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD);
440:         DMGetLocalSection(dmErrAux, &sec);
441:         PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD);
442:       }
443:       DMViewFromOptions(dmErrAux, NULL, "-dm_err_view");
444:       ISViewFromOptions(subis[0], NULL, "-super_is_view");
445:       ISViewFromOptions(subis[1], NULL, "-super_is_view");
446:       DMGetGlobalVector(dmErrAux, &uErr);
447:       VecViewFromOptions(u, NULL, "-map_vec_view");
448:       VecViewFromOptions(uAdjProj, NULL, "-map_vec_view");
449:       VecViewFromOptions(uErr, NULL, "-map_vec_view");
450:       VecISCopy(uErr, subis[0], SCATTER_FORWARD, u);
451:       VecISCopy(uErr, subis[1], SCATTER_FORWARD, uAdjProj);
452:       DMRestoreGlobalVector(dm, &uAdjProj);
453:       for (i = 0; i < 2; ++i) ISDestroy(&subis[i]);
454:       PetscFree(subis);
455:       DMGetLocalVector(dmErrAux, &uErrLoc);
456:       DMGlobalToLocalBegin(dm, uErr, INSERT_VALUES, uErrLoc);
457:       DMGlobalToLocalEnd(dm, uErr, INSERT_VALUES, uErrLoc);
458:       DMRestoreGlobalVector(dmErrAux, &uErr);
459:       DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, uErrLoc);
460:       /*   Compute cellwise error estimate */
461:       VecSet(errorEst, 0.0);
462:       DMPlexComputeCellwiseIntegralFEM(dmAdj, uAdj, errorEst, &user);
463:       DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, NULL);
464:       DMRestoreLocalVector(dmErrAux, &uErrLoc);
465:       DMDestroy(&dmErrAux);
466:       /*   Plot cellwise error vector */
467:       VecViewFromOptions(errorEst, NULL, "-error_view");
468:       /*   Compute ratio of estimate (sum over cells) with actual L_2 error */
469:       DMComputeL2Diff(dm, 0.0, funcs, ctxs, u, &errorL2Norm);
470:       DMPlexComputeL2DiffVec(dm, 0.0, funcs, ctxs, u, errorL2);
471:       VecViewFromOptions(errorL2, NULL, "-l2_error_view");
472:       VecNorm(errorL2, NORM_INFINITY, &errorL2Tot);
473:       VecNorm(errorEst, NORM_INFINITY, &errorEstTot);
474:       VecGetSize(errorEst, &N);
475:       VecPointwiseDivide(errorEst, errorEst, errorL2);
476:       PetscObjectSetName((PetscObject)errorEst, "Error ratio");
477:       VecViewFromOptions(errorEst, NULL, "-error_ratio_view");
478:       PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " L2 error: %g Error Ratio: %g/%g = %g\n", N, (double)errorL2Norm, (double)errorEstTot, (double)PetscSqrtReal(errorL2Tot), (double)(errorEstTot / PetscSqrtReal(errorL2Tot)));
479:       DMRestoreGlobalVector(dmErr, &errorEst);
480:       DMRestoreGlobalVector(dmErr, &errorL2);
481:       DMDestroy(&dmErr);
482:     }
483:     DMDestroy(&dmAdj);
484:     VecDestroy(&uAdj);
485:     SNESDestroy(&snesAdj);
486:   }
487:   /* Cleanup */
488:   VecDestroy(&u);
489:   SNESDestroy(&snes);
490:   DMDestroy(&dm);
491:   PetscFinalize();
492:   return 0;
493: }

495: /*TEST

497:   test:
498:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
499:     suffix: 2d_p1_conv
500:     requires: triangle
501:     args: -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
502:   test:
503:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
504:     suffix: 2d_p2_conv
505:     requires: triangle
506:     args: -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
507:   test:
508:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
509:     suffix: 2d_p3_conv
510:     requires: triangle
511:     args: -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
512:   test:
513:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
514:     suffix: 2d_q1_conv
515:     args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
516:   test:
517:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
518:     suffix: 2d_q2_conv
519:     args: -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
520:   test:
521:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
522:     suffix: 2d_q3_conv
523:     args: -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
524:   test:
525:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
526:     suffix: 2d_q1_shear_conv
527:     args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
528:   test:
529:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
530:     suffix: 2d_q2_shear_conv
531:     args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
532:   test:
533:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
534:     suffix: 2d_q3_shear_conv
535:     args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
536:   test:
537:     # Using -convest_num_refine 3 we get L_2 convergence rate: 1.7
538:     suffix: 3d_p1_conv
539:     requires: ctetgen
540:     args: -dm_plex_dim 3 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1
541:   test:
542:     # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 2.8
543:     suffix: 3d_p2_conv
544:     requires: ctetgen
545:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1
546:   test:
547:     # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 4.0
548:     suffix: 3d_p3_conv
549:     requires: ctetgen
550:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1
551:   test:
552:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.8
553:     suffix: 3d_q1_conv
554:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1
555:   test:
556:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.8
557:     suffix: 3d_q2_conv
558:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1
559:   test:
560:     # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 3.8
561:     suffix: 3d_q3_conv
562:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1
563:   test:
564:     suffix: 2d_p1_fas_full
565:     requires: triangle
566:     args: -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \
567:       -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full -snes_fas_full_total \
568:         -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \
569:         -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \
570:           -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \
571:             -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10
572:   test:
573:     suffix: 2d_p1_fas_full_homogeneous
574:     requires: triangle
575:     args: -homogeneous -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \
576:       -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full \
577:         -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \
578:         -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \
579:           -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \
580:             -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10

582:   test:
583:     suffix: 2d_p1_scalable
584:     requires: triangle
585:     args: -potential_petscspace_degree 1 -dm_refine 3 \
586:       -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned \
587:       -pc_type gamg -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 \
588:         -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \
589:         -pc_gamg_coarse_eq_limit 1000 \
590:         -pc_gamg_threshold 0.05 \
591:         -pc_gamg_threshold_scale .0 \
592:         -mg_levels_ksp_type chebyshev \
593:         -mg_levels_ksp_max_it 1 \
594:         -mg_levels_pc_type jacobi \
595:       -matptap_via scalable
596:   test:
597:     suffix: 2d_p1_gmg_vcycle
598:     requires: triangle
599:     args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
600:           -ksp_rtol 5e-10 -pc_type mg \
601:             -mg_levels_ksp_max_it 1 \
602:             -mg_levels_esteig_ksp_type cg \
603:             -mg_levels_esteig_ksp_max_it 10 \
604:             -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
605:             -mg_levels_pc_type jacobi
606:   test:
607:     suffix: 2d_p1_gmg_fcycle
608:     requires: triangle
609:     args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
610:           -ksp_rtol 5e-10 -pc_type mg -pc_mg_type full \
611:             -mg_levels_ksp_max_it 2 \
612:             -mg_levels_esteig_ksp_type cg \
613:             -mg_levels_esteig_ksp_max_it 10 \
614:             -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
615:             -mg_levels_pc_type jacobi
616:   test:
617:     suffix: 2d_p1_gmg_vcycle_adapt
618:     requires: triangle
619:     args: -petscpartitioner_type simple -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
620:           -ksp_rtol 5e-10 -pc_type mg -pc_mg_galerkin -pc_mg_adapt_interp_coarse_space harmonic -pc_mg_adapt_interp_n 8 \
621:             -mg_levels_ksp_max_it 1 \
622:             -mg_levels_esteig_ksp_type cg \
623:             -mg_levels_esteig_ksp_max_it 10 \
624:             -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
625:             -mg_levels_pc_type jacobi
626:   test:
627:     suffix: 2d_p1_spectral_0
628:     requires: triangle fftw !complex
629:     args: -dm_plex_box_faces 1,1 -potential_petscspace_degree 1 -dm_refine 6 -spectral -fft_view
630:   test:
631:     suffix: 2d_p1_spectral_1
632:     requires: triangle fftw !complex
633:     nsize: 2
634:     args: -dm_plex_box_faces 4,4 -potential_petscspace_degree 1 -spectral -fft_view
635:   test:
636:     suffix: 2d_p1_adj_0
637:     requires: triangle
638:     args: -potential_petscspace_degree 1 -dm_refine 1 -adjoint -adjoint_petscspace_degree 1 -error_petscspace_degree 0
639:   test:
640:     nsize: 2
641:     requires: kokkos_kernels
642:     suffix: kokkos
643:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,3,6 -petscpartitioner_type simple -dm_plex_simplex 0 -potential_petscspace_degree 1 \
644:          -dm_refine 0 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -pc_type gamg -pc_gamg_coarse_eq_limit 1000 -pc_gamg_threshold 0.0 \
645:          -pc_gamg_threshold_scale .5 -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 2 -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 \
646:          -ksp_monitor -snes_monitor -dm_view -dm_mat_type aijkokkos -dm_vec_type kokkos

648: TEST*/