Actual source code: ex36.c

      1: static char help[] = "First example in homogenization book\n\n";

      3: #include <petscsnes.h>
      4: #include <petscdmplex.h>
      5: #include <petscds.h>
      6: #include <petscconvest.h>
      7: #include <petscbag.h>

      9: /*
     10:   To control the refinement use -dm_plex_box_faces <n> or -dm_refine <k>, or both

     12:   To see the exact and computed solutions

     14:     -compare_view draw -draw_size 500,500 -draw_pause -1

     16:   To see the delay in convergence of the discretization use

     18:     -snes_convergence_estimate -convest_num_refine 7 -convest_monitor

     20:   and to see the proper rate use

     22:     -dm_refine 5 -snes_convergence_estimate -convest_num_refine 2 -convest_monitor
     23: */

     25: typedef enum {
     26:   MOD_CONSTANT,
     27:   MOD_OSCILLATORY,
     28:   MOD_TANH,
     29:   NUM_MOD_TYPES
     30: } ModType;
     31: const char *modTypes[NUM_MOD_TYPES + 1] = {"constant", "oscillatory", "tanh", "unknown"};

     33: /* Constants */
     34: enum {
     35:   EPSILON,
     36:   NUM_CONSTANTS
     37: };

     39: typedef struct {
     40:   PetscReal epsilon; /* Wavelength of fine scale oscillation */
     41: } Parameter;

     43: typedef struct {
     44:   PetscBag bag;     /* Holds problem parameters */
     45:   ModType  modType; /* Model type */
     46: } AppCtx;

     48: static PetscErrorCode trig_homogeneous_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
     49: {
     50:   PetscInt d;
     51:   *u = 1.0;
     52:   for (d = 0; d < dim; ++d) *u *= PetscSinReal(2.0 * PETSC_PI * x[d]);
     53:   return PETSC_SUCCESS;
     54: }

     56: static PetscErrorCode oscillatory_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
     57: {
     58:   Parameter      *param = (Parameter *)ctx;
     59:   const PetscReal eps   = param->epsilon;

     61:   u[0] = x[0] - x[0] * x[0] + (eps / (2. * PETSC_PI)) * (0.5 - x[0]) * PetscSinReal(2. * PETSC_PI * x[0] / eps) + PetscSqr(eps / (2. * PETSC_PI)) * (1. - PetscCosReal(2. * PETSC_PI * x[0] / eps));
     62:   return PETSC_SUCCESS;
     63: }

     65: 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[])
     66: {
     67:   PetscInt d;
     68:   for (d = 0; d < dim; ++d) {
     69:     PetscScalar v = 1.;
     70:     for (PetscInt e = 0; e < dim; e++) {
     71:       if (e == d) {
     72:         v *= -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
     73:       } else {
     74:         v *= PetscSinReal(2.0 * PETSC_PI * x[d]);
     75:       }
     76:     }
     77:     f0[0] += v;
     78:   }
     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 void f0_oscillatory_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[])
     94: {
     95:   f0[0] = -1.;
     96: }

     98: static void f1_oscillatory_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[])
     99: {
    100:   const PetscReal eps = PetscRealPart(constants[EPSILON]);

    102:   f1[0] = u_x[0] / (2. + PetscCosReal(2. * PETSC_PI * x[0] / eps));
    103: }

    105: static void g3_oscillatory_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[])
    106: {
    107:   const PetscReal eps = PetscRealPart(constants[EPSILON]);

    109:   g3[0] = 1. / (2. + PetscCosReal(2. * PETSC_PI * x[0] / eps));
    110: }

    112: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
    113: {
    114:   PetscInt mod;

    116:   PetscFunctionBeginUser;
    117:   options->modType = MOD_CONSTANT;
    118:   PetscOptionsBegin(comm, "", "Homogenization Problem Options", "DMPLEX");
    119:   mod = options->modType;
    120:   PetscCall(PetscOptionsEList("-mod_type", "The model type", "ex36.c", modTypes, NUM_MOD_TYPES, modTypes[options->modType], &mod, NULL));
    121:   options->modType = (ModType)mod;
    122:   PetscOptionsEnd();
    123:   PetscFunctionReturn(PETSC_SUCCESS);
    124: }

    126: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *user)
    127: {
    128:   PetscBag   bag;
    129:   Parameter *p;

    131:   PetscFunctionBeginUser;
    132:   PetscCall(PetscBagCreate(comm, sizeof(Parameter), &user->bag));
    133:   PetscCall(PetscBagGetData(user->bag, &p));
    134:   PetscCall(PetscBagSetName(user->bag, "par", "Homogenization parameters"));
    135:   bag = user->bag;
    136:   PetscCall(PetscBagRegisterReal(bag, &p->epsilon, 1.0, "epsilon", "Wavelength of fine scale oscillation"));
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
    141: {
    142:   PetscFunctionBeginUser;
    143:   PetscCall(DMCreate(comm, dm));
    144:   PetscCall(DMSetType(*dm, DMPLEX));
    145:   PetscCall(DMSetFromOptions(*dm));
    146:   PetscCall(DMSetApplicationContext(*dm, user));
    147:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    148:   PetscFunctionReturn(PETSC_SUCCESS);
    149: }

    151: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
    152: {
    153:   PetscDS             ds;
    154:   DMLabel             label;
    155:   PetscSimplePointFn *ex;
    156:   const PetscInt      id = 1;
    157:   void               *ctx;

    159:   PetscFunctionBeginUser;
    160:   PetscCall(DMGetDS(dm, &ds));
    161:   PetscCall(PetscBagGetData(user->bag, &ctx));
    162:   switch (user->modType) {
    163:   case MOD_CONSTANT:
    164:     PetscCall(PetscDSSetResidual(ds, 0, f0_trig_homogeneous_u, f1_u));
    165:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
    166:     PetscCall(DMGetLabel(dm, "marker", &label));
    167:     ex = trig_homogeneous_u;
    168:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)ex, NULL, ctx, NULL));
    169:     break;
    170:   case MOD_OSCILLATORY:
    171:     PetscCall(PetscDSSetResidual(ds, 0, f0_oscillatory_u, f1_oscillatory_u));
    172:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_oscillatory_uu));
    173:     PetscCall(DMGetLabel(dm, "marker", &label));
    174:     ex = oscillatory_u;
    175:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)ex, NULL, ctx, NULL));
    176:     break;
    177:   default:
    178:     SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported model type: %s (%d)", modTypes[PetscMin(user->modType, NUM_MOD_TYPES)], user->modType);
    179:   }
    180:   PetscCall(PetscDSSetExactSolution(ds, 0, ex, ctx));
    181:   /* Setup constants */
    182:   {
    183:     Parameter  *param;
    184:     PetscScalar constants[NUM_CONSTANTS];

    186:     PetscCall(PetscBagGetData(user->bag, &param));

    188:     constants[EPSILON] = param->epsilon;
    189:     PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
    190:   }
    191:   PetscFunctionReturn(PETSC_SUCCESS);
    192: }

    194: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
    195: {
    196:   DM        cdm = dm;
    197:   PetscFE   fe;
    198:   PetscBool simplex;
    199:   PetscInt  dim;
    200:   char      prefix[PETSC_MAX_PATH_LEN];

    202:   PetscFunctionBeginUser;
    203:   PetscCall(DMGetDimension(dm, &dim));
    204:   PetscCall(DMPlexIsSimplex(dm, &simplex));
    205:   /* Create finite element */
    206:   PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
    207:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe));
    208:   PetscCall(PetscObjectSetName((PetscObject)fe, name));
    209:   /* Set discretization and boundary conditions for each mesh */
    210:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
    211:   PetscCall(DMCreateDS(dm));
    212:   PetscCall((*setup)(dm, user));
    213:   while (cdm) {
    214:     PetscCall(DMCopyDisc(dm, cdm));
    215:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    216:   }
    217:   PetscCall(PetscFEDestroy(&fe));
    218:   PetscFunctionReturn(PETSC_SUCCESS);
    219: }

    221: static PetscErrorCode CompareView(Vec u)
    222: {
    223:   DM                dm;
    224:   Vec               v[2], lv[2], exact;
    225:   PetscOptions      options;
    226:   PetscViewer       viewer;
    227:   PetscViewerFormat format;
    228:   PetscBool         flg;
    229:   PetscInt          i;
    230:   const char       *name, *prefix;

    232:   PetscFunctionBeginUser;
    233:   PetscCall(VecGetDM(u, &dm));
    234:   PetscCall(PetscObjectGetOptions((PetscObject)dm, &options));
    235:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    236:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)dm), options, prefix, "-compare_view", &viewer, &format, &flg));
    237:   if (flg) {
    238:     PetscCall(DMGetGlobalVector(dm, &exact));
    239:     PetscCall(DMComputeExactSolution(dm, 0.0, exact, NULL));
    240:     v[0] = u;
    241:     v[1] = exact;
    242:     for (i = 0; i < 2; ++i) {
    243:       PetscCall(DMGetLocalVector(dm, &lv[i]));
    244:       PetscCall(PetscObjectGetName((PetscObject)v[i], &name));
    245:       PetscCall(PetscObjectSetName((PetscObject)lv[i], name));
    246:       PetscCall(DMGlobalToLocalBegin(dm, v[i], INSERT_VALUES, lv[i]));
    247:       PetscCall(DMGlobalToLocalEnd(dm, v[i], INSERT_VALUES, lv[i]));
    248:       PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, lv[i], 0., NULL, NULL, NULL));
    249:     }
    250:     PetscCall(DMPlexVecView1D(dm, 2, lv, viewer));
    251:     for (i = 0; i < 2; ++i) PetscCall(DMRestoreLocalVector(dm, &lv[i]));
    252:     PetscCall(DMRestoreGlobalVector(dm, &exact));
    253:     PetscCall(PetscViewerDestroy(&viewer));
    254:   }
    255:   PetscFunctionReturn(PETSC_SUCCESS);
    256: }

    258: typedef struct {
    259:   Mat Mcoarse;   /* Mass matrix on the coarse space */
    260:   Mat Mfine;     /* Mass matrix on the fine space */
    261:   Mat Ifine;     /* Interpolator from coarse to fine */
    262:   Vec Iscale;    /* Scaling vector for restriction */
    263:   KSP kspCoarse; /* Solver for the coarse mass matrix */
    264:   Vec tmpfine;   /* Temporary vector in the fine space */
    265:   Vec tmpcoarse; /* Temporary vector in the coarse space */
    266: } ProjStruct;

    268: static PetscErrorCode DestroyCoarseProjection(Mat Pi)
    269: {
    270:   ProjStruct *ctx;

    272:   PetscFunctionBegin;
    273:   PetscCall(MatShellGetContext(Pi, &ctx));
    274:   PetscCall(MatDestroy(&ctx->Mcoarse));
    275:   PetscCall(MatDestroy(&ctx->Mfine));
    276:   PetscCall(MatDestroy(&ctx->Ifine));
    277:   PetscCall(VecDestroy(&ctx->Iscale));
    278:   PetscCall(KSPDestroy(&ctx->kspCoarse));
    279:   PetscCall(VecDestroy(&ctx->tmpcoarse));
    280:   PetscCall(VecDestroy(&ctx->tmpfine));
    281:   PetscCall(PetscFree(ctx));
    282:   PetscCall(MatShellSetContext(Pi, NULL));
    283:   PetscFunctionReturn(PETSC_SUCCESS);
    284: }

    286: static PetscErrorCode CoarseProjection(Mat Pi, Vec x, Vec y)
    287: {
    288:   ProjStruct *ctx;

    290:   PetscFunctionBegin;
    291:   PetscCall(MatShellGetContext(Pi, &ctx));
    292:   PetscCall(MatMult(ctx->Mfine, x, ctx->tmpfine));
    293:   PetscCall(PetscObjectSetName((PetscObject)ctx->tmpfine, "Fine DG RHS"));
    294:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->tmpfine, "fine_dg_"));
    295:   PetscCall(VecViewFromOptions(ctx->tmpfine, NULL, "-rhs_view"));
    296:   PetscCall(MatMultTranspose(ctx->Ifine, ctx->tmpfine, ctx->tmpcoarse));
    297:   PetscCall(PetscObjectSetName((PetscObject)ctx->tmpcoarse, "Coarse DG RHS"));
    298:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->tmpcoarse, "coarse_dg_"));
    299:   PetscCall(VecViewFromOptions(ctx->tmpcoarse, NULL, "-rhs_view"));
    300:   PetscCall(VecPointwiseMult(ctx->tmpcoarse, ctx->Iscale, ctx->tmpcoarse));
    301:   PetscCall(VecViewFromOptions(ctx->tmpcoarse, NULL, "-rhs_view"));
    302:   PetscCall(KSPSolve(ctx->kspCoarse, ctx->tmpcoarse, y));
    303:   PetscFunctionReturn(PETSC_SUCCESS);
    304: }

    306: static PetscErrorCode CreateCoarseProjection(DM dmc, DM dmf, Mat *Pi)
    307: {
    308:   ProjStruct *ctx;
    309:   PetscInt    m, n, M, N;

    311:   PetscFunctionBegin;
    312:   PetscCall(PetscMalloc1(1, &ctx));
    313:   PetscCall(DMCreateGlobalVector(dmc, &ctx->tmpcoarse));
    314:   PetscCall(DMCreateGlobalVector(dmf, &ctx->tmpfine));
    315:   PetscCall(VecGetLocalSize(ctx->tmpcoarse, &m));
    316:   PetscCall(VecGetSize(ctx->tmpcoarse, &M));
    317:   PetscCall(VecGetLocalSize(ctx->tmpfine, &n));
    318:   PetscCall(VecGetSize(ctx->tmpfine, &N));
    319:   PetscCall(DMCreateMassMatrix(dmc, dmc, &ctx->Mcoarse));
    320:   PetscCall(DMCreateMassMatrix(dmf, dmf, &ctx->Mfine));
    321:   PetscCall(DMCreateInterpolation(dmc, dmf, &ctx->Ifine, &ctx->Iscale));
    322:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)dmc), &ctx->kspCoarse));
    323:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->kspCoarse, "coarse_"));
    324:   PetscCall(KSPSetOperators(ctx->kspCoarse, ctx->Mcoarse, ctx->Mcoarse));
    325:   PetscCall(KSPSetFromOptions(ctx->kspCoarse));
    326:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)dmc), m, n, M, N, ctx, Pi));
    327:   PetscCall(MatShellSetOperation(*Pi, MATOP_DESTROY, (PetscErrorCodeFn *)DestroyCoarseProjection));
    328:   PetscCall(MatShellSetOperation(*Pi, MATOP_MULT, (PetscErrorCodeFn *)CoarseProjection));
    329:   PetscFunctionReturn(PETSC_SUCCESS);
    330: }

    332: typedef struct {
    333:   Mat Ifdg; /* Embed the fine space into the DG version */
    334:   Mat Pi;   /* The L_2 stable projection to the DG coarse space */
    335:   Vec tmpc; /* A temporary vector in the DG coarse space */
    336:   Vec tmpf; /* A temporary vector in the DG fine space */
    337: } QuasiInterp;

    339: static PetscErrorCode DestroyQuasiInterpolator(Mat P)
    340: {
    341:   QuasiInterp *ctx;

    343:   PetscFunctionBegin;
    344:   PetscCall(MatShellGetContext(P, &ctx));
    345:   PetscCall(MatDestroy(&ctx->Ifdg));
    346:   PetscCall(MatDestroy(&ctx->Pi));
    347:   PetscCall(VecDestroy(&ctx->tmpc));
    348:   PetscCall(VecDestroy(&ctx->tmpf));
    349:   PetscCall(PetscFree(ctx));
    350:   PetscCall(MatShellSetContext(P, NULL));
    351:   PetscFunctionReturn(PETSC_SUCCESS);
    352: }

    354: static PetscErrorCode QuasiInterpolate(Mat P, Vec x, Vec y)
    355: {
    356:   QuasiInterp *ctx;
    357:   DM           dmcdg, dmc;
    358:   Vec          ly;

    360:   PetscFunctionBegin;
    361:   PetscCall(MatShellGetContext(P, &ctx));
    362:   PetscCall(MatMult(ctx->Ifdg, x, ctx->tmpf));

    364:   PetscCall(PetscObjectSetName((PetscObject)ctx->tmpf, "Fine DG Potential"));
    365:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->tmpf, "fine_dg_"));
    366:   PetscCall(VecViewFromOptions(ctx->tmpf, NULL, "-vec_view"));
    367:   PetscCall(MatMult(ctx->Pi, x, ctx->tmpc));

    369:   PetscCall(PetscObjectSetName((PetscObject)ctx->tmpc, "Coarse DG Potential"));
    370:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->tmpc, "coarse_dg_"));
    371:   PetscCall(VecViewFromOptions(ctx->tmpc, NULL, "-vec_view"));
    372:   PetscCall(VecGetDM(ctx->tmpc, &dmcdg));

    374:   PetscCall(VecGetDM(y, &dmc));
    375:   PetscCall(DMGetLocalVector(dmc, &ly));
    376:   PetscCall(DMPlexComputeClementInterpolant(dmcdg, ctx->tmpc, ly));
    377:   PetscCall(DMLocalToGlobalBegin(dmc, ly, INSERT_VALUES, y));
    378:   PetscCall(DMLocalToGlobalEnd(dmc, ly, INSERT_VALUES, y));
    379:   PetscCall(DMRestoreLocalVector(dmc, &ly));
    380:   PetscFunctionReturn(PETSC_SUCCESS);
    381: }

    383: static PetscErrorCode CreateQuasiInterpolator(DM dmc, DM dmf, Mat *P)
    384: {
    385:   QuasiInterp   *ctx;
    386:   DM             dmcdg, dmfdg;
    387:   PetscFE        fe;
    388:   Vec            x, y;
    389:   DMPolytopeType ct;
    390:   PetscInt       dim, cStart, m, n, M, N;

    392:   PetscFunctionBegin;
    393:   PetscCall(PetscCalloc1(1, &ctx));
    394:   PetscCall(DMGetGlobalVector(dmc, &x));
    395:   PetscCall(DMGetGlobalVector(dmf, &y));
    396:   PetscCall(VecGetLocalSize(x, &m));
    397:   PetscCall(VecGetSize(x, &M));
    398:   PetscCall(VecGetLocalSize(y, &n));
    399:   PetscCall(VecGetSize(y, &N));
    400:   PetscCall(DMRestoreGlobalVector(dmc, &x));
    401:   PetscCall(DMRestoreGlobalVector(dmf, &y));

    403:   PetscCall(DMClone(dmf, &dmfdg));
    404:   PetscCall(DMGetDimension(dmfdg, &dim));
    405:   PetscCall(DMPlexGetHeightStratum(dmfdg, 0, &cStart, NULL));
    406:   PetscCall(DMPlexGetCellType(dmfdg, cStart, &ct));
    407:   PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, "fine_dg_", PETSC_DETERMINE, &fe));
    408:   PetscCall(DMSetField(dmfdg, 0, NULL, (PetscObject)fe));
    409:   PetscCall(PetscFEDestroy(&fe));
    410:   PetscCall(DMCreateDS(dmfdg));
    411:   PetscCall(DMCreateInterpolation(dmf, dmfdg, &ctx->Ifdg, NULL));
    412:   PetscCall(DMCreateGlobalVector(dmfdg, &ctx->tmpf));
    413:   PetscCall(DMDestroy(&dmfdg));

    415:   PetscCall(DMClone(dmc, &dmcdg));
    416:   PetscCall(DMGetDimension(dmcdg, &dim));
    417:   PetscCall(DMPlexGetHeightStratum(dmcdg, 0, &cStart, NULL));
    418:   PetscCall(DMPlexGetCellType(dmcdg, cStart, &ct));
    419:   PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, "coarse_dg_", PETSC_DETERMINE, &fe));
    420:   PetscCall(DMSetField(dmcdg, 0, NULL, (PetscObject)fe));
    421:   PetscCall(PetscFEDestroy(&fe));
    422:   PetscCall(DMCreateDS(dmcdg));

    424:   PetscCall(CreateCoarseProjection(dmcdg, dmf, &ctx->Pi));
    425:   PetscCall(DMCreateGlobalVector(dmcdg, &ctx->tmpc));
    426:   PetscCall(DMDestroy(&dmcdg));

    428:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)dmc), m, n, M, N, ctx, P));
    429:   PetscCall(MatShellSetOperation(*P, MATOP_DESTROY, (PetscErrorCodeFn *)DestroyQuasiInterpolator));
    430:   PetscCall(MatShellSetOperation(*P, MATOP_MULT, (PetscErrorCodeFn *)QuasiInterpolate));
    431:   PetscFunctionReturn(PETSC_SUCCESS);
    432: }

    434: static PetscErrorCode CoarseTest(DM dm, Vec u, AppCtx *user)
    435: {
    436:   DM  dmc;
    437:   Mat P; /* The quasi-interpolator to the coarse space */
    438:   Vec uc;

    440:   PetscFunctionBegin;
    441:   if (user->modType == MOD_CONSTANT) PetscFunctionReturn(PETSC_SUCCESS);
    442:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &dmc));
    443:   PetscCall(DMSetType(dmc, DMPLEX));
    444:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dmc, "coarse_"));
    445:   PetscCall(DMSetApplicationContext(dmc, user));
    446:   PetscCall(DMSetFromOptions(dmc));
    447:   PetscCall(DMViewFromOptions(dmc, NULL, "-dm_view"));

    449:   PetscCall(SetupDiscretization(dmc, "potential", SetupPrimalProblem, user));
    450:   PetscCall(DMCreateGlobalVector(dmc, &uc));
    451:   PetscCall(PetscObjectSetName((PetscObject)uc, "potential"));
    452:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)uc, "coarse_"));

    454:   PetscCall(CreateQuasiInterpolator(dmc, dm, &P));
    455: #if 1
    456:   PetscCall(MatMult(P, u, uc));
    457: #else
    458:   {
    459:     Mat In;
    460:     Vec sc;

    462:     PetscCall(DMCreateInterpolation(dmc, dm, &In, &sc));
    463:     PetscCall(MatMultTranspose(In, u, uc));
    464:     PetscCall(VecPointwiseMult(uc, sc, uc));
    465:     PetscCall(MatDestroy(&In));
    466:     PetscCall(VecDestroy(&sc));
    467:   }
    468: #endif
    469:   PetscCall(CompareView(uc));

    471:   PetscCall(MatDestroy(&P));
    472:   PetscCall(VecDestroy(&uc));
    473:   PetscCall(DMDestroy(&dmc));
    474:   PetscFunctionReturn(PETSC_SUCCESS);
    475: }

    477: int main(int argc, char **argv)
    478: {
    479:   DM     dm;   /* Problem specification */
    480:   SNES   snes; /* Nonlinear solver */
    481:   Vec    u;    /* Solutions */
    482:   AppCtx user; /* User-defined work context */

    484:   PetscFunctionBeginUser;
    485:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    486:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    487:   PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
    488:   /* Primal system */
    489:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    490:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    491:   PetscCall(SNESSetDM(snes, dm));
    492:   PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
    493:   PetscCall(DMCreateGlobalVector(dm, &u));
    494:   PetscCall(VecSet(u, 0.0));
    495:   PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
    496:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
    497:   PetscCall(SNESSetFromOptions(snes));
    498:   PetscCall(DMSNESCheckFromOptions(snes, u));
    499:   PetscCall(SNESSolve(snes, NULL, u));
    500:   PetscCall(SNESGetSolution(snes, &u));
    501:   PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
    502:   PetscCall(CompareView(u));
    503:   /* Looking at a coarse problem */
    504:   PetscCall(CoarseTest(dm, u, &user));
    505:   /* Cleanup */
    506:   PetscCall(VecDestroy(&u));
    507:   PetscCall(SNESDestroy(&snes));
    508:   PetscCall(DMDestroy(&dm));
    509:   PetscCall(PetscBagDestroy(&user.bag));
    510:   PetscCall(PetscFinalize());
    511:   return 0;
    512: }

    514: /*TEST

    516:   test:
    517:     suffix: 1d_p1_constant
    518:     args: -dm_plex_dim 1 -dm_plex_box_faces 4 -potential_petscspace_degree 1 -dmsnes_check

    520:   test:
    521:     suffix: 1d_p1_constant_conv
    522:     args: -dm_plex_dim 1 -dm_plex_box_faces 4 -potential_petscspace_degree 1 \
    523:           -snes_convergence_estimate -convest_num_refine 2

    525:   test:
    526:     suffix: 1d_p1_oscillatory
    527:     args: -mod_type oscillatory -epsilon 0.03125 \
    528:           -dm_plex_dim 1 -dm_plex_box_faces 4 -potential_petscspace_degree 1 -dm_refine 2 -dmsnes_check \
    529:           -coarse_dm_plex_dim 1 -coarse_dm_plex_box_faces 4 -coarse_dm_plex_hash_location \
    530:           -fine_dg_petscspace_degree 1 -fine_dg_petscdualspace_lagrange_continuity 0 \
    531:           -coarse_dg_petscspace_degree 1 -coarse_dg_petscdualspace_lagrange_continuity 0

    533: TEST*/