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 <petscdmceed.h>
  9: #include <petscsnes.h>
 10: #include <petscds.h>
 11: #include <petscconvest.h>

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

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

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

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

 44: /* Compute integral of (residual of solution)*(adjoint solution - projection of adjoint solution) */
 45: 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[])
 46: {
 47:   obj[0] = a[aOff[0]] * (u[0] - a[aOff[1]]);
 48: }

 50: 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[])
 51: {
 52:   PetscInt d;
 53:   for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
 54: }

 56: 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[])
 57: {
 58:   PetscInt d;
 59:   for (d = 0; d < dim; ++d) {
 60:     PetscScalar v = 1.;
 61:     for (PetscInt e = 0; e < dim; e++) {
 62:       if (e == d) {
 63:         v *= -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
 64:       } else {
 65:         v *= PetscSinReal(2.0 * PETSC_PI * x[d]);
 66:       }
 67:     }
 68:     f0[0] += v;
 69:   }
 70: }

 72: 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[])
 73: {
 74:   f0[0] = 1.0;
 75: }

 77: 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[])
 78: {
 79:   f0[0] = a[0];
 80: }

 82: 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[])
 83: {
 84:   PetscInt d;
 85:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
 86: }

 88: 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[])
 89: {
 90:   PetscInt d;
 91:   for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
 92: }

 94: PLEXFE_QFUNCTION(Laplace, f0_trig_inhomogeneous_u, f1_u)

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

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

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

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

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

139:       PetscCall(PetscSectionGetOffset(coordSection, v, &off));
140:       if (PetscAbsReal(planeCoord[p] - PetscRealPart(coords[off + planeDir[p]])) < PETSC_SMALL) PetscCall(DMLabelSetValue(label, v, 1));
141:     }
142:   }
143:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
144:   PetscFunctionReturn(PETSC_SUCCESS);
145: }

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

160:     PetscCall(CreateSpectralPlanes(*dm, 2, planeDir, planeCoord, user));
161:   }
162:   PetscFunctionReturn(PETSC_SUCCESS);
163: }

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

173:   PetscFunctionBeginUser;
174:   PetscCall(DMGetDS(dm, &ds));
175:   PetscCall(PetscDSSetResidual(ds, 0, f0, f1_u));
176:   PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
177:   PetscCall(PetscDSSetExactSolution(ds, 0, ex, user));
178:   PetscCall(DMGetLabel(dm, "marker", &label));
179:   if (label) PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))ex, NULL, user, NULL));
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

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

189:   PetscFunctionBeginUser;
190:   PetscCall(DMGetDS(dm, &ds));
191:   PetscCall(PetscDSSetResidual(ds, 0, f0_unity_u, f1_u));
192:   PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
193:   PetscCall(PetscDSSetObjective(ds, 0, obj_error_u));
194:   PetscCall(DMGetLabel(dm, "marker", &label));
195:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))zero, NULL, user, NULL));
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: static PetscErrorCode SetupErrorProblem(DM dm, AppCtx *user)
200: {
201:   PetscDS prob;

203:   PetscFunctionBeginUser;
204:   PetscCall(DMGetDS(dm, &prob));
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

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

217:   PetscFunctionBeginUser;
218:   PetscCall(DMGetDimension(dm, &dim));
219:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
220:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
221:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
222:   /* Create finite element */
223:   PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
224:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe));
225:   PetscCall(PetscObjectSetName((PetscObject)fe, name));
226:   /* Set discretization and boundary conditions for each mesh */
227:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
228:   PetscCall(DMCreateDS(dm));
229:   PetscCall((*setup)(dm, user));
230:   while (cdm) {
231:     PetscCall(DMCopyDisc(dm, cdm));
232:     /* TODO: Check whether the boundary of coarse meshes is marked */
233:     PetscCall(DMGetCoarseDM(cdm, &cdm));
234:   }
235:   PetscCall(PetscFEDestroy(&fe));
236: #ifdef PETSC_HAVE_LIBCEED
237:   PetscBool useCeed;
238:   PetscCall(DMPlexGetUseCeed(dm, &useCeed));
239:   if (useCeed) PetscCall(DMCeedCreate(dm, PETSC_TRUE, PlexQFunctionLaplace, PlexQFunctionLaplace_loc));
240: #endif
241:   PetscFunctionReturn(PETSC_SUCCESS);
242: }

244: static PetscErrorCode ComputeSpectral(Vec u, PetscInt numPlanes, const PetscInt planeDir[], const PetscReal planeCoord[], AppCtx *user)
245: {
246:   MPI_Comm           comm;
247:   DM                 dm;
248:   PetscSection       coordSection, section;
249:   Vec                coordinates, uloc;
250:   const PetscScalar *coords, *array;
251:   PetscInt           p;
252:   PetscMPIInt        size, rank;

254:   PetscFunctionBeginUser;
255:   if (!user->spectral) PetscFunctionReturn(PETSC_SUCCESS);
256:   PetscCall(VecGetDM(u, &dm));
257:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
258:   PetscCallMPI(MPI_Comm_size(comm, &size));
259:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
260:   PetscCall(DMGetLocalVector(dm, &uloc));
261:   PetscCall(DMGlobalToLocalBegin(dm, u, INSERT_VALUES, uloc));
262:   PetscCall(DMGlobalToLocalEnd(dm, u, INSERT_VALUES, uloc));
263:   PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, uloc, 0.0, NULL, NULL, NULL));
264:   PetscCall(VecViewFromOptions(uloc, NULL, "-sol_view"));
265:   PetscCall(DMGetLocalSection(dm, &section));
266:   PetscCall(VecGetArrayRead(uloc, &array));
267:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
268:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
269:   PetscCall(VecGetArrayRead(coordinates, &coords));
270:   for (p = 0; p < numPlanes; ++p) {
271:     DMLabel         label;
272:     char            name[PETSC_MAX_PATH_LEN];
273:     Mat             F;
274:     Vec             x, y;
275:     IS              stratum;
276:     PetscReal      *ray, *gray;
277:     PetscScalar    *rvals, *svals, *gsvals;
278:     PetscInt       *perm, *nperm;
279:     PetscInt        n, N, i, j, off, offu;
280:     const PetscInt *points;

282:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p));
283:     PetscCall(DMGetLabel(dm, name, &label));
284:     PetscCall(DMLabelGetStratumIS(label, 1, &stratum));
285:     PetscCall(ISGetLocalSize(stratum, &n));
286:     PetscCall(ISGetIndices(stratum, &points));
287:     PetscCall(PetscMalloc2(n, &ray, n, &svals));
288:     for (i = 0; i < n; ++i) {
289:       PetscCall(PetscSectionGetOffset(coordSection, points[i], &off));
290:       PetscCall(PetscSectionGetOffset(section, points[i], &offu));
291:       ray[i]   = PetscRealPart(coords[off + ((planeDir[p] + 1) % 2)]);
292:       svals[i] = array[offu];
293:     }
294:     /* Gather the ray data to proc 0 */
295:     if (size > 1) {
296:       PetscMPIInt *cnt, *displs, p;

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

354: static PetscErrorCode ComputeAdjoint(Vec u, AppCtx *user)
355: {
356:   PetscFunctionBegin;
357:   if (!user->adjoint) PetscFunctionReturn(PETSC_SUCCESS);
358:   DM   dm, dmAdj;
359:   SNES snesAdj;
360:   Vec  uAdj;

362:   PetscCall(VecGetDM(u, &dm));
363:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snesAdj));
364:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)snesAdj, "adjoint_"));
365:   PetscCall(DMClone(dm, &dmAdj));
366:   PetscCall(SNESSetDM(snesAdj, dmAdj));
367:   PetscCall(SetupDiscretization(dmAdj, "adjoint", SetupAdjointProblem, user));
368:   PetscCall(DMCreateGlobalVector(dmAdj, &uAdj));
369:   PetscCall(VecSet(uAdj, 0.0));
370:   PetscCall(PetscObjectSetName((PetscObject)uAdj, "adjoint"));
371:   PetscCall(DMPlexSetSNESLocalFEM(dmAdj, PETSC_FALSE, &user));
372:   PetscCall(SNESSetFromOptions(snesAdj));
373:   PetscCall(SNESSolve(snesAdj, NULL, uAdj));
374:   PetscCall(SNESGetSolution(snesAdj, &uAdj));
375:   PetscCall(VecViewFromOptions(uAdj, NULL, "-adjoint_view"));
376:   /* Error representation */
377:   {
378:     DM        dmErr, dmErrAux, dms[2];
379:     Vec       errorEst, errorL2, uErr, uErrLoc, uAdjLoc, uAdjProj;
380:     IS       *subis;
381:     PetscReal errorEstTot, errorL2Norm, errorL2Tot;
382:     PetscInt  N, i;
383:     PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {user->homogeneous ? trig_homogeneous_u : trig_inhomogeneous_u};
384:     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};
385:     void *ctxs[1] = {0};

387:     ctxs[0] = user;
388:     PetscCall(DMClone(dm, &dmErr));
389:     PetscCall(SetupDiscretization(dmErr, "error", SetupErrorProblem, user));
390:     PetscCall(DMGetGlobalVector(dmErr, &errorEst));
391:     PetscCall(DMGetGlobalVector(dmErr, &errorL2));
392:     /*   Compute auxiliary data (solution and projection of adjoint solution) */
393:     PetscCall(DMGetLocalVector(dmAdj, &uAdjLoc));
394:     PetscCall(DMGlobalToLocalBegin(dmAdj, uAdj, INSERT_VALUES, uAdjLoc));
395:     PetscCall(DMGlobalToLocalEnd(dmAdj, uAdj, INSERT_VALUES, uAdjLoc));
396:     PetscCall(DMGetGlobalVector(dm, &uAdjProj));
397:     PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, uAdjLoc));
398:     PetscCall(DMProjectField(dm, 0.0, u, identity, INSERT_VALUES, uAdjProj));
399:     PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL));
400:     PetscCall(DMRestoreLocalVector(dmAdj, &uAdjLoc));
401:     /*   Attach auxiliary data */
402:     dms[0] = dm;
403:     dms[1] = dm;
404:     PetscCall(DMCreateSuperDM(dms, 2, &subis, &dmErrAux));
405:     if (0) {
406:       PetscSection sec;

408:       PetscCall(DMGetLocalSection(dms[0], &sec));
409:       PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD));
410:       PetscCall(DMGetLocalSection(dms[1], &sec));
411:       PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD));
412:       PetscCall(DMGetLocalSection(dmErrAux, &sec));
413:       PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD));
414:     }
415:     PetscCall(DMViewFromOptions(dmErrAux, NULL, "-dm_err_view"));
416:     PetscCall(ISViewFromOptions(subis[0], NULL, "-super_is_view"));
417:     PetscCall(ISViewFromOptions(subis[1], NULL, "-super_is_view"));
418:     PetscCall(DMGetGlobalVector(dmErrAux, &uErr));
419:     PetscCall(VecViewFromOptions(u, NULL, "-map_vec_view"));
420:     PetscCall(VecViewFromOptions(uAdjProj, NULL, "-map_vec_view"));
421:     PetscCall(VecViewFromOptions(uErr, NULL, "-map_vec_view"));
422:     PetscCall(VecISCopy(uErr, subis[0], SCATTER_FORWARD, u));
423:     PetscCall(VecISCopy(uErr, subis[1], SCATTER_FORWARD, uAdjProj));
424:     PetscCall(DMRestoreGlobalVector(dm, &uAdjProj));
425:     for (i = 0; i < 2; ++i) PetscCall(ISDestroy(&subis[i]));
426:     PetscCall(PetscFree(subis));
427:     PetscCall(DMGetLocalVector(dmErrAux, &uErrLoc));
428:     PetscCall(DMGlobalToLocalBegin(dm, uErr, INSERT_VALUES, uErrLoc));
429:     PetscCall(DMGlobalToLocalEnd(dm, uErr, INSERT_VALUES, uErrLoc));
430:     PetscCall(DMRestoreGlobalVector(dmErrAux, &uErr));
431:     PetscCall(DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, uErrLoc));
432:     /*   Compute cellwise error estimate */
433:     PetscCall(VecSet(errorEst, 0.0));
434:     PetscCall(DMPlexComputeCellwiseIntegralFEM(dmAdj, uAdj, errorEst, user));
435:     PetscCall(DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, NULL));
436:     PetscCall(DMRestoreLocalVector(dmErrAux, &uErrLoc));
437:     PetscCall(DMDestroy(&dmErrAux));
438:     /*   Plot cellwise error vector */
439:     PetscCall(VecViewFromOptions(errorEst, NULL, "-error_view"));
440:     /*   Compute ratio of estimate (sum over cells) with actual L_2 error */
441:     PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, u, &errorL2Norm));
442:     PetscCall(DMPlexComputeL2DiffVec(dm, 0.0, funcs, ctxs, u, errorL2));
443:     PetscCall(VecViewFromOptions(errorL2, NULL, "-l2_error_view"));
444:     PetscCall(VecNorm(errorL2, NORM_INFINITY, &errorL2Tot));
445:     PetscCall(VecNorm(errorEst, NORM_INFINITY, &errorEstTot));
446:     PetscCall(VecGetSize(errorEst, &N));
447:     PetscCall(VecPointwiseDivide(errorEst, errorEst, errorL2));
448:     PetscCall(PetscObjectSetName((PetscObject)errorEst, "Error ratio"));
449:     PetscCall(VecViewFromOptions(errorEst, NULL, "-error_ratio_view"));
450:     PetscCall(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))));
451:     PetscCall(DMRestoreGlobalVector(dmErr, &errorEst));
452:     PetscCall(DMRestoreGlobalVector(dmErr, &errorL2));
453:     PetscCall(DMDestroy(&dmErr));
454:   }
455:   PetscCall(DMDestroy(&dmAdj));
456:   PetscCall(VecDestroy(&uAdj));
457:   PetscCall(SNESDestroy(&snesAdj));
458:   PetscFunctionReturn(PETSC_SUCCESS);
459: }

461: static PetscErrorCode ErrorView(Vec u, AppCtx *user)
462: {
463:   PetscErrorCode (*sol)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
464:   void     *ctx;
465:   DM        dm;
466:   PetscDS   ds;
467:   PetscReal error;
468:   PetscInt  N;

470:   PetscFunctionBegin;
471:   if (!user->viewError) PetscFunctionReturn(PETSC_SUCCESS);
472:   PetscCall(VecGetDM(u, &dm));
473:   PetscCall(DMGetDS(dm, &ds));
474:   PetscCall(PetscDSGetExactSolution(ds, 0, &sol, &ctx));
475:   PetscCall(VecGetSize(u, &N));
476:   PetscCall(DMComputeL2Diff(dm, 0.0, &sol, &ctx, u, &error));
477:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " L2 error: %g\n", N, (double)error));
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: int main(int argc, char **argv)
482: {
483:   DM        dm;   /* Problem specification */
484:   SNES      snes; /* Nonlinear solver */
485:   Vec       u;    /* Solutions */
486:   AppCtx    user; /* User-defined work context */
487:   PetscInt  planeDir[2]   = {0, 1};
488:   PetscReal planeCoord[2] = {0., 1.};

490:   PetscFunctionBeginUser;
491:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
492:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
493:   /* Primal system */
494:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
495:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
496:   PetscCall(SNESSetDM(snes, dm));
497:   PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
498:   PetscCall(DMCreateGlobalVector(dm, &u));
499:   PetscCall(VecSet(u, 0.0));
500:   PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
501:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
502:   PetscCall(SNESSetFromOptions(snes));
503:   PetscCall(SNESSolve(snes, NULL, u));
504:   PetscCall(SNESGetSolution(snes, &u));
505:   PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
506:   PetscCall(ErrorView(u, &user));
507:   PetscCall(ComputeSpectral(u, 2, planeDir, planeCoord, &user));
508:   PetscCall(ComputeAdjoint(u, &user));
509:   /* Cleanup */
510:   PetscCall(VecDestroy(&u));
511:   PetscCall(SNESDestroy(&snes));
512:   PetscCall(DMDestroy(&dm));
513:   PetscCall(PetscFinalize());
514:   return 0;
515: }

517: /*TEST

519:   test:
520:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
521:     suffix: 2d_p1_conv
522:     requires: triangle
523:     args: -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
524:   test:
525:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
526:     suffix: 2d_p2_conv
527:     requires: triangle
528:     args: -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
529:   test:
530:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
531:     suffix: 2d_p3_conv
532:     requires: triangle
533:     args: -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
534:   test:
535:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
536:     suffix: 2d_q1_conv
537:     args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
538:   test:
539:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
540:     suffix: 2d_q2_conv
541:     args: -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
542:   test:
543:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
544:     suffix: 2d_q3_conv
545:     args: -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
546:   test:
547:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
548:     suffix: 2d_q1_ceed_conv
549:     requires: libceed
550:     args: -dm_plex_use_ceed -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
551:   test:
552:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
553:     suffix: 2d_q2_ceed_conv
554:     requires: libceed
555:     args: -dm_plex_use_ceed -dm_plex_simplex 0 -potential_petscspace_degree 2 -cdm_default_quadrature_order 2 \
556:           -snes_convergence_estimate -convest_num_refine 2
557:   test:
558:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
559:     suffix: 2d_q3_ceed_conv
560:     requires: libceed
561:     args: -dm_plex_use_ceed -dm_plex_simplex 0 -potential_petscspace_degree 3 -cdm_default_quadrature_order 3 \
562:           -snes_convergence_estimate -convest_num_refine 2
563:   test:
564:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
565:     suffix: 2d_q1_shear_conv
566:     args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
567:   test:
568:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
569:     suffix: 2d_q2_shear_conv
570:     args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
571:   test:
572:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
573:     suffix: 2d_q3_shear_conv
574:     args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
575:   test:
576:     # Using -convest_num_refine 3 we get L_2 convergence rate: 1.7
577:     suffix: 3d_p1_conv
578:     requires: ctetgen
579:     args: -dm_plex_dim 3 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1
580:   test:
581:     # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 2.8
582:     suffix: 3d_p2_conv
583:     requires: ctetgen
584:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1
585:   test:
586:     # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 4.0
587:     suffix: 3d_p3_conv
588:     requires: ctetgen
589:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1
590:   test:
591:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.8
592:     suffix: 3d_q1_conv
593:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1
594:   test:
595:     # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.8
596:     suffix: 3d_q2_conv
597:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1
598:   test:
599:     # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 3.8
600:     suffix: 3d_q3_conv
601:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1
602:   test:
603:     suffix: 2d_p1_fas_full
604:     requires: triangle
605:     args: -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \
606:       -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full -snes_fas_full_total \
607:         -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \
608:         -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \
609:           -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \
610:             -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10
611:   test:
612:     suffix: 2d_p1_fas_full_homogeneous
613:     requires: triangle
614:     args: -homogeneous -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \
615:       -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full \
616:         -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \
617:         -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \
618:           -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \
619:             -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10

621:   test:
622:     suffix: 2d_p1_scalable
623:     requires: triangle
624:     args: -potential_petscspace_degree 1 -dm_refine 3 \
625:       -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned \
626:       -pc_type gamg -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 \
627:         -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \
628:         -pc_gamg_coarse_eq_limit 1000 \
629:         -pc_gamg_threshold 0.05 \
630:         -pc_gamg_threshold_scale .0 \
631:         -mg_levels_ksp_type chebyshev \
632:         -mg_levels_ksp_max_it 1 \
633:         -mg_levels_pc_type jacobi \
634:       -matptap_via scalable
635:   test:
636:     suffix: 2d_p1_gmg_vcycle
637:     requires: triangle
638:     args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
639:           -ksp_rtol 5e-10 -pc_type mg \
640:             -mg_levels_ksp_max_it 1 \
641:             -mg_levels_esteig_ksp_type cg \
642:             -mg_levels_esteig_ksp_max_it 10 \
643:             -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
644:             -mg_levels_pc_type jacobi
645:   testset:
646:     args: -potential_petscspace_degree 1 -dm_refine_hierarchy 3 \
647:           -ksp_rtol 5e-10 -pc_type mg -pc_mg_type full \
648:             -mg_levels_ksp_max_it 2 \
649:             -mg_levels_esteig_ksp_type cg \
650:             -mg_levels_esteig_ksp_max_it 10 \
651:             -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
652:             -mg_levels_pc_type jacobi
653:     test:
654:       suffix: 2d_p1_gmg_fcycle
655:       requires: triangle
656:       args: -dm_plex_box_faces 2,2
657:     test:
658:       suffix: 2d_q1_gmg_fcycle
659:       args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2
660:     test:
661:       suffix: 3d_p1_gmg_fcycle
662:       requires: ctetgen
663:       args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2
664:     test:
665:       suffix: 3d_q1_gmg_fcycle
666:       args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2
667:   test:
668:     suffix: 2d_p1_gmg_vcycle_adapt
669:     requires: triangle
670:     args: -petscpartitioner_type simple -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
671:           -ksp_rtol 5e-10 -pc_type mg -pc_mg_galerkin -pc_mg_adapt_interp_coarse_space harmonic -pc_mg_adapt_interp_n 8 \
672:             -mg_levels_ksp_max_it 1 \
673:             -mg_levels_esteig_ksp_type cg \
674:             -mg_levels_esteig_ksp_max_it 10 \
675:             -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
676:             -mg_levels_pc_type jacobi
677:   test:
678:     suffix: 2d_p1_spectral_0
679:     requires: triangle fftw !complex
680:     args: -dm_plex_box_faces 1,1 -potential_petscspace_degree 1 -dm_refine 6 -spectral -fft_view
681:   test:
682:     suffix: 2d_p1_spectral_1
683:     requires: triangle fftw !complex
684:     nsize: 2
685:     args: -dm_plex_box_faces 4,4 -potential_petscspace_degree 1 -spectral -fft_view
686:   test:
687:     suffix: 2d_p1_adj_0
688:     requires: triangle
689:     args: -potential_petscspace_degree 1 -dm_refine 1 -adjoint -adjoint_petscspace_degree 1 -error_petscspace_degree 0
690:   test:
691:     nsize: 2
692:     requires: kokkos_kernels
693:     suffix: kokkos
694:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,3,6 -petscpartitioner_type simple -dm_plex_simplex 0 -potential_petscspace_degree 1 \
695:          -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 \
696:          -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 \
697:          -ksp_monitor -snes_monitor -dm_view -dm_mat_type aijkokkos -dm_vec_type kokkos

699: TEST*/