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 boundary 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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 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: {
 95:   PetscFunctionBeginUser;
 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:   PetscCall(PetscOptionsBool("-shear", "Shear the domain", "ex13.c", options->shear, &options->shear, NULL));
104:   PetscCall(PetscOptionsBool("-spectral", "Look at the spectrum along planes of the solution", "ex13.c", options->spectral, &options->spectral, NULL));
105:   PetscCall(PetscOptionsBool("-adjoint", "Solve the adjoint problem", "ex13.c", options->adjoint, &options->adjoint, NULL));
106:   PetscCall(PetscOptionsBool("-homogeneous", "Use homogeneous boundary conditions", "ex13.c", options->homogeneous, &options->homogeneous, NULL));
107:   PetscCall(PetscOptionsBool("-error_view", "Output the solution error", "ex13.c", options->viewError, &options->viewError, NULL));
108:   PetscOptionsEnd();
109:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

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

157:     PetscCall(CreateSpectralPlanes(*dm, 2, planeDir, planeCoord, user));
158:   }
159:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

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

200:   PetscFunctionBeginUser;
201:   PetscCall(DMGetDS(dm, &prob));
202:   PetscFunctionReturn(PETSC_SUCCESS);
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];

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

245:   PetscFunctionBeginUser;
246:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
247:   PetscCallMPI(MPI_Comm_size(comm, &size));
248:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
249:   PetscCall(DMGetLocalVector(dm, &uloc));
250:   PetscCall(DMGlobalToLocalBegin(dm, u, INSERT_VALUES, uloc));
251:   PetscCall(DMGlobalToLocalEnd(dm, u, INSERT_VALUES, uloc));
252:   PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, uloc, 0.0, NULL, NULL, NULL));
253:   PetscCall(VecViewFromOptions(uloc, NULL, "-sol_view"));
254:   PetscCall(DMGetLocalSection(dm, &section));
255:   PetscCall(VecGetArrayRead(uloc, &array));
256:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
257:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
258:   PetscCall(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:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p));
272:     PetscCall(DMGetLabel(dm, name, &label));
273:     PetscCall(DMLabelGetStratumIS(label, 1, &stratum));
274:     PetscCall(ISGetLocalSize(stratum, &n));
275:     PetscCall(ISGetIndices(stratum, &points));
276:     PetscCall(PetscMalloc2(n, &ray, n, &svals));
277:     for (i = 0; i < n; ++i) {
278:       PetscCall(PetscSectionGetOffset(coordSection, points[i], &off));
279:       PetscCall(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:       PetscCall(PetscCalloc2(size, &cnt, size, &displs));
288:       PetscCallMPI(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:       PetscCall(PetscMalloc2(N, &gray, N, &gsvals));
292:       PetscCallMPI(MPI_Gatherv(ray, n, MPIU_REAL, gray, cnt, displs, MPIU_REAL, 0, comm));
293:       PetscCallMPI(MPI_Gatherv(svals, n, MPIU_SCALAR, gsvals, cnt, displs, MPIU_SCALAR, 0, comm));
294:       PetscCall(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:       PetscCall(PetscMalloc2(N, &perm, N, &nperm));
303:       for (i = 0; i < N; ++i) perm[i] = i;
304:       PetscCall(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:       PetscCall(MatCreateFFT(PETSC_COMM_SELF, 1, &j, MATFFTW, &F));
312:       PetscCall(MatCreateVecs(F, &x, &y));
313:       PetscCall(PetscObjectSetName((PetscObject)y, name));
314:       PetscCall(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:       PetscCall(PetscFree2(perm, nperm));
321:       if (size > 1) PetscCall(PetscFree2(gray, gsvals));
322:       PetscCall(VecRestoreArray(x, &rvals));
323:       /* Do FFT along the ray */
324:       PetscCall(MatMult(F, x, y));
325:       /* Chop FFT */
326:       PetscCall(VecChop(y, PETSC_SMALL));
327:       PetscCall(VecViewFromOptions(x, NULL, "-real_view"));
328:       PetscCall(VecViewFromOptions(y, NULL, "-fft_view"));
329:       PetscCall(VecDestroy(&x));
330:       PetscCall(VecDestroy(&y));
331:       PetscCall(MatDestroy(&F));
332:     }
333:     PetscCall(ISRestoreIndices(stratum, &points));
334:     PetscCall(ISDestroy(&stratum));
335:     PetscCall(PetscFree2(ray, svals));
336:   }
337:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
338:   PetscCall(VecRestoreArrayRead(uloc, &array));
339:   PetscCall(DMRestoreLocalVector(dm, &uloc));
340:   PetscFunctionReturn(PETSC_SUCCESS);
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 */

350:   PetscFunctionBeginUser;
351:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
352:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
353:   /* Primal system */
354:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
355:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
356:   PetscCall(SNESSetDM(snes, dm));
357:   PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
358:   PetscCall(DMCreateGlobalVector(dm, &u));
359:   PetscCall(VecSet(u, 0.0));
360:   PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
361:   PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user));
362:   PetscCall(SNESSetFromOptions(snes));
363:   PetscCall(SNESSolve(snes, NULL, u));
364:   PetscCall(SNESGetSolution(snes, &u));
365:   PetscCall(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:     PetscCall(DMGetDS(dm, &ds));
374:     PetscCall(PetscDSGetExactSolution(ds, 0, &sol, &ctx));
375:     PetscCall(VecGetSize(u, &N));
376:     PetscCall(DMComputeL2Diff(dm, 0.0, &sol, &ctx, u, &error));
377:     PetscCall(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:     PetscCall(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:     PetscCall(SNESCreate(PETSC_COMM_WORLD, &snesAdj));
392:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)snesAdj, "adjoint_"));
393:     PetscCall(DMClone(dm, &dmAdj));
394:     PetscCall(SNESSetDM(snesAdj, dmAdj));
395:     PetscCall(SetupDiscretization(dmAdj, "adjoint", SetupAdjointProblem, &user));
396:     PetscCall(DMCreateGlobalVector(dmAdj, &uAdj));
397:     PetscCall(VecSet(uAdj, 0.0));
398:     PetscCall(PetscObjectSetName((PetscObject)uAdj, "adjoint"));
399:     PetscCall(DMPlexSetSNESLocalFEM(dmAdj, &user, &user, &user));
400:     PetscCall(SNESSetFromOptions(snesAdj));
401:     PetscCall(SNESSolve(snesAdj, NULL, uAdj));
402:     PetscCall(SNESGetSolution(snesAdj, &uAdj));
403:     PetscCall(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:       PetscCall(DMClone(dm, &dmErr));
417:       PetscCall(SetupDiscretization(dmErr, "error", SetupErrorProblem, &user));
418:       PetscCall(DMGetGlobalVector(dmErr, &errorEst));
419:       PetscCall(DMGetGlobalVector(dmErr, &errorL2));
420:       /*   Compute auxiliary data (solution and projection of adjoint solution) */
421:       PetscCall(DMGetLocalVector(dmAdj, &uAdjLoc));
422:       PetscCall(DMGlobalToLocalBegin(dmAdj, uAdj, INSERT_VALUES, uAdjLoc));
423:       PetscCall(DMGlobalToLocalEnd(dmAdj, uAdj, INSERT_VALUES, uAdjLoc));
424:       PetscCall(DMGetGlobalVector(dm, &uAdjProj));
425:       PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, uAdjLoc));
426:       PetscCall(DMProjectField(dm, 0.0, u, identity, INSERT_VALUES, uAdjProj));
427:       PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL));
428:       PetscCall(DMRestoreLocalVector(dmAdj, &uAdjLoc));
429:       /*   Attach auxiliary data */
430:       dms[0] = dm;
431:       dms[1] = dm;
432:       PetscCall(DMCreateSuperDM(dms, 2, &subis, &dmErrAux));
433:       if (0) {
434:         PetscSection sec;

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