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:     PetscMPIInt     in;
    281:     const PetscInt *points;

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

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

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

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

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

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

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

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

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

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

    519: /*TEST

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

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

    702: TEST*/