Actual source code: ex45.c

      1: static char help[] = "Heat Equation in 2d and 3d with finite elements.\n\
      2: We solve the heat equation in a rectangular\n\
      3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
      4: Contributed by: Julian Andrej <juan@tf.uni-kiel.de>\n\n\n";

      6: #include <petscdmplex.h>
      7: #include <petscds.h>
      8: #include <petscts.h>

     10: /*
     11:   Heat equation:

     13:     du/dt - \Delta u + f = 0
     14: */

     16: typedef enum {
     17:   SOL_QUADRATIC_LINEAR,
     18:   SOL_QUADRATIC_TRIG,
     19:   SOL_TRIG_LINEAR,
     20:   SOL_TRIG_TRIG,
     21:   NUM_SOLUTION_TYPES
     22: } SolutionType;
     23: const char *solutionTypes[NUM_SOLUTION_TYPES + 1] = {"quadratic_linear", "quadratic_trig", "trig_linear", "trig_trig", "unknown"};

     25: typedef struct {
     26:   SolutionType solType; /* Type of exact solution */
     27:   /* Solver setup */
     28:   PetscBool expTS;        /* Flag for explicit timestepping */
     29:   PetscBool lumped;       /* Lump the mass matrix */
     30:   PetscInt  remesh_every; /* Remesh every number of steps */
     31:   DM        remesh_dm;    /* New DM after remeshing */
     32: } AppCtx;

     34: /*
     35: Exact 2D solution:
     36:   u    = 2t + x^2 + y^2
     37:   u_t  = 2
     38:   \Delta u = 2 + 2 = 4
     39:   f    = 2
     40:   F(u) = 2 - (2 + 2) + 2 = 0

     42: Exact 3D solution:
     43:   u = 3t + x^2 + y^2 + z^2
     44:   F(u) = 3 - (2 + 2 + 2) + 3 = 0
     45: */
     46: static PetscErrorCode mms_quad_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
     47: {
     48:   PetscInt d;

     50:   *u = dim * time;
     51:   for (d = 0; d < dim; ++d) *u += x[d] * x[d];
     52:   return PETSC_SUCCESS;
     53: }

     55: static PetscErrorCode mms_quad_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
     56: {
     57:   *u = dim;
     58:   return PETSC_SUCCESS;
     59: }

     61: static void f0_quad_lin_exp(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[])
     62: {
     63:   f0[0] = -(PetscScalar)dim;
     64: }
     65: static void f0_quad_lin(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:   PetscScalar exp[1] = {0.};
     68:   f0_quad_lin_exp(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, exp);
     69:   f0[0] = u_t[0] - exp[0];
     70: }

     72: /*
     73: Exact 2D solution:
     74:   u = 2*cos(t) + x^2 + y^2
     75:   F(u) = -2*sint(t) - (2 + 2) + 2*sin(t) + 4 = 0

     77: Exact 3D solution:
     78:   u = 3*cos(t) + x^2 + y^2 + z^2
     79:   F(u) = -3*sin(t) - (2 + 2 + 2) + 3*sin(t) + 6 = 0
     80: */
     81: static PetscErrorCode mms_quad_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
     82: {
     83:   PetscInt d;

     85:   *u = dim * PetscCosReal(time);
     86:   for (d = 0; d < dim; ++d) *u += x[d] * x[d];
     87:   return PETSC_SUCCESS;
     88: }

     90: static PetscErrorCode mms_quad_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
     91: {
     92:   *u = -dim * PetscSinReal(time);
     93:   return PETSC_SUCCESS;
     94: }

     96: static void f0_quad_trig_exp(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[])
     97: {
     98:   f0[0] = -dim * (PetscSinReal(t) + 2.0);
     99: }
    100: static void f0_quad_trig(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[])
    101: {
    102:   PetscScalar exp[1] = {0.};
    103:   f0_quad_trig_exp(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, exp);
    104:   f0[0] = u_t[0] - exp[0];
    105: }

    107: /*
    108: Exact 2D solution:
    109:   u = 2\pi^2 t + cos(\pi x) + cos(\pi y)
    110:   F(u) = 2\pi^2 - \pi^2 (cos(\pi x) + cos(\pi y)) + \pi^2 (cos(\pi x) + cos(\pi y)) - 2\pi^2 = 0

    112: Exact 3D solution:
    113:   u = 3\pi^2 t + cos(\pi x) + cos(\pi y) + cos(\pi z)
    114:   F(u) = 3\pi^2 - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) + \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) - 3\pi^2 = 0
    115: */
    116: static PetscErrorCode mms_trig_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
    117: {
    118:   PetscInt d;

    120:   *u = dim * PetscSqr(PETSC_PI) * time;
    121:   for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI * x[d]);
    122:   return PETSC_SUCCESS;
    123: }

    125: static PetscErrorCode mms_trig_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
    126: {
    127:   *u = dim * PetscSqr(PETSC_PI);
    128:   return PETSC_SUCCESS;
    129: }

    131: static void f0_trig_lin(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[])
    132: {
    133:   PetscInt d;
    134:   f0[0] = u_t[0];
    135:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI) * (PetscCosReal(PETSC_PI * x[d]) - 1.0);
    136: }

    138: /*
    139: Exact 2D solution:
    140:   u    = pi^2 cos(t) + cos(\pi x) + cos(\pi y)
    141:   u_t  = -pi^2 sin(t)
    142:   \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y))
    143:   f    = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y))
    144:   F(u) = -\pi^2 sin(t) + \pi^2 (cos(\pi x) + cos(\pi y)) - \pi^2 (cos(\pi x) + cos(\pi y)) + \pi^2 sin(t) = 0

    146: Exact 3D solution:
    147:   u    = pi^2 cos(t) + cos(\pi x) + cos(\pi y) + cos(\pi z)
    148:   u_t  = -pi^2 sin(t)
    149:   \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z))
    150:   f    = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z))
    151:   F(u) = -\pi^2 sin(t) + \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) + \pi^2 sin(t) = 0
    152: */
    153: static PetscErrorCode mms_trig_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
    154: {
    155:   PetscInt d;

    157:   *u = PetscSqr(PETSC_PI) * PetscCosReal(time);
    158:   for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI * x[d]);
    159:   return PETSC_SUCCESS;
    160: }

    162: static PetscErrorCode mms_trig_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
    163: {
    164:   *u = -PetscSqr(PETSC_PI) * PetscSinReal(time);
    165:   return PETSC_SUCCESS;
    166: }

    168: static void f0_trig_trig_exp(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[])
    169: {
    170:   PetscInt d;
    171:   f0[0] -= PetscSqr(PETSC_PI) * PetscSinReal(t);
    172:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI) * PetscCosReal(PETSC_PI * x[d]);
    173: }
    174: static void f0_trig_trig(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[])
    175: {
    176:   PetscScalar exp[1] = {0.};
    177:   f0_trig_trig_exp(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, exp);
    178:   f0[0] = u_t[0] - exp[0];
    179: }

    181: static void f1_temp_exp(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[])
    182: {
    183:   for (PetscInt d = 0; d < dim; ++d) f1[d] = -u_x[d];
    184: }
    185: static void f1_temp(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[])
    186: {
    187:   for (PetscInt d = 0; d < dim; ++d) f1[d] = u_x[d];
    188: }

    190: static void g3_temp(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[])
    191: {
    192:   for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
    193: }

    195: static void g0_temp(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 g0[])
    196: {
    197:   g0[0] = u_tShift * 1.0;
    198: }

    200: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
    201: {
    202:   PetscInt sol;

    204:   PetscFunctionBeginUser;
    205:   options->solType      = SOL_QUADRATIC_LINEAR;
    206:   options->expTS        = PETSC_FALSE;
    207:   options->lumped       = PETSC_TRUE;
    208:   options->remesh_every = 0;

    210:   PetscOptionsBegin(comm, "", "Heat Equation Options", "DMPLEX");
    211:   sol = options->solType;
    212:   PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex45.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL));
    213:   options->solType = (SolutionType)sol;
    214:   PetscCall(PetscOptionsBool("-explicit", "Use explicit timestepping", "ex45.c", options->expTS, &options->expTS, NULL));
    215:   PetscCall(PetscOptionsBool("-lumped", "Lump the mass matrix", "ex45.c", options->lumped, &options->lumped, NULL));
    216:   PetscCall(PetscOptionsInt("-remesh_every", "Remesh every number of steps", "ex45.c", options->remesh_every, &options->remesh_every, NULL));
    217:   PetscOptionsEnd();
    218:   PetscFunctionReturn(PETSC_SUCCESS);
    219: }

    221: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
    222: {
    223:   PetscFunctionBeginUser;
    224:   PetscCall(DMCreate(comm, dm));
    225:   PetscCall(DMSetType(*dm, DMPLEX));
    226:   PetscCall(DMSetFromOptions(*dm));
    227:   {
    228:     char      convType[256];
    229:     PetscBool flg;
    230:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
    231:     PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", __FILE__, DMList, DMPLEX, convType, 256, &flg));
    232:     PetscOptionsEnd();
    233:     if (flg) {
    234:       DM dmConv;
    235:       PetscCall(DMConvert(*dm, convType, &dmConv));
    236:       if (dmConv) {
    237:         PetscCall(DMDestroy(dm));
    238:         *dm = dmConv;
    239:         PetscCall(DMSetFromOptions(*dm));
    240:         PetscCall(DMSetUp(*dm));
    241:       }
    242:     }
    243:   }
    244:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    245:   PetscFunctionReturn(PETSC_SUCCESS);
    246: }

    248: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
    249: {
    250:   PetscDS        ds;
    251:   DMLabel        label;
    252:   const PetscInt id = 1;

    254:   PetscFunctionBeginUser;
    255:   PetscCall(DMGetLabel(dm, "marker", &label));
    256:   PetscCall(DMGetDS(dm, &ds));
    257:   PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_temp, NULL, NULL, g3_temp));
    258:   switch (ctx->solType) {
    259:   case SOL_QUADRATIC_LINEAR:
    260:     PetscCall(PetscDSSetResidual(ds, 0, f0_quad_lin, f1_temp));
    261:     PetscCall(PetscDSSetRHSResidual(ds, 0, f0_quad_lin_exp, f1_temp_exp));
    262:     PetscCall(PetscDSSetExactSolution(ds, 0, mms_quad_lin, ctx));
    263:     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_lin_t, ctx));
    264:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)mms_quad_lin, (PetscVoidFn *)mms_quad_lin_t, ctx, NULL));
    265:     break;
    266:   case SOL_QUADRATIC_TRIG:
    267:     PetscCall(PetscDSSetResidual(ds, 0, f0_quad_trig, f1_temp));
    268:     PetscCall(PetscDSSetRHSResidual(ds, 0, f0_quad_trig_exp, f1_temp_exp));
    269:     PetscCall(PetscDSSetExactSolution(ds, 0, mms_quad_trig, ctx));
    270:     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_trig_t, ctx));
    271:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)mms_quad_trig, (PetscVoidFn *)mms_quad_trig_t, ctx, NULL));
    272:     break;
    273:   case SOL_TRIG_LINEAR:
    274:     PetscCall(PetscDSSetResidual(ds, 0, f0_trig_lin, f1_temp));
    275:     PetscCall(PetscDSSetExactSolution(ds, 0, mms_trig_lin, ctx));
    276:     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_lin_t, ctx));
    277:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)mms_trig_lin, (PetscVoidFn *)mms_trig_lin_t, ctx, NULL));
    278:     break;
    279:   case SOL_TRIG_TRIG:
    280:     PetscCall(PetscDSSetResidual(ds, 0, f0_trig_trig, f1_temp));
    281:     PetscCall(PetscDSSetRHSResidual(ds, 0, f0_trig_trig_exp, f1_temp_exp));
    282:     PetscCall(PetscDSSetExactSolution(ds, 0, mms_trig_trig, ctx));
    283:     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_trig_t, ctx));
    284:     break;
    285:   default:
    286:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType);
    287:   }
    288:   PetscFunctionReturn(PETSC_SUCCESS);
    289: }

    291: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
    292: {
    293:   DM        plex, cdm = dm;
    294:   PetscFE   fe;
    295:   PetscBool simplex;
    296:   PetscInt  dim;

    298:   PetscFunctionBeginUser;
    299:   PetscCall(DMGetDimension(dm, &dim));
    300:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    301:   PetscCall(DMPlexIsSimplex(plex, &simplex));
    302:   PetscCall(DMDestroy(&plex));
    303:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", -1, &fe));
    304:   PetscCall(PetscObjectSetName((PetscObject)fe, "temperature"));
    305:   /* Set discretization and boundary conditions for each mesh */
    306:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
    307:   PetscCall(DMCreateDS(dm));
    308:   if (ctx->expTS) {
    309:     PetscDS ds;

    311:     PetscCall(DMGetDS(dm, &ds));
    312:     PetscCall(PetscDSSetImplicit(ds, 0, PETSC_FALSE));
    313:   }
    314:   PetscCall(SetupProblem(dm, ctx));
    315:   while (cdm) {
    316:     PetscCall(DMCopyDisc(dm, cdm));
    317:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    318:   }
    319:   PetscCall(PetscFEDestroy(&fe));
    320:   PetscFunctionReturn(PETSC_SUCCESS);
    321: }

    323: #include <petsc/private/dmpleximpl.h>
    324: static PetscErrorCode Remesh(DM dm, Vec U, DM *newdm)
    325: {
    326:   PetscFunctionBeginUser;
    327:   PetscCall(DMViewFromOptions(dm, NULL, "-remesh_dmin_view"));
    328:   *newdm = NULL;

    330:   PetscInt  dim;
    331:   DM        plex;
    332:   PetscBool simplex;
    333:   PetscCall(DMGetDimension(dm, &dim));
    334:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    335:   PetscCall(DMPlexIsSimplex(plex, &simplex));
    336:   PetscCall(DMDestroy(&plex));

    338:   DM      dmGrad;
    339:   PetscFE fe;
    340:   PetscCall(DMClone(dm, &dmGrad));
    341:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "grad_", -1, &fe));
    342:   PetscCall(DMSetField(dmGrad, 0, NULL, (PetscObject)fe));
    343:   PetscCall(PetscFEDestroy(&fe));
    344:   PetscCall(DMCreateDS(dmGrad));

    346:   Vec locU, locG;
    347:   PetscCall(DMGetLocalVector(dm, &locU));
    348:   PetscCall(DMGetLocalVector(dmGrad, &locG));
    349:   PetscCall(DMGlobalToLocal(dm, U, INSERT_VALUES, locU));
    350:   PetscCall(VecViewFromOptions(locU, NULL, "-remesh_input_view"));
    351:   PetscCall(DMPlexComputeGradientClementInterpolant(dm, locU, locG));
    352:   PetscCall(VecViewFromOptions(locG, NULL, "-remesh_grad_view"));

    354:   const PetscScalar *g;
    355:   PetscScalar       *u;
    356:   PetscInt           n;
    357:   PetscCall(VecGetLocalSize(locU, &n));
    358:   PetscCall(VecGetArrayWrite(locU, &u));
    359:   PetscCall(VecGetArrayRead(locG, &g));
    360:   for (PetscInt i = 0; i < n; i++) {
    361:     PetscReal norm = 0.0;
    362:     for (PetscInt d = 0; d < dim; d++) norm += PetscSqr(PetscRealPart(g[dim * i + d]));
    363:     u[i] = PetscSqrtReal(norm);
    364:   }
    365:   PetscCall(VecRestoreArrayRead(locG, &g));
    366:   PetscCall(VecRestoreArrayWrite(locU, &u));

    368:   DM  dmM;
    369:   Vec metric;
    370:   PetscCall(DMClone(dm, &dmM));
    371:   PetscCall(DMPlexMetricCreateIsotropic(dmM, 0, locU, &metric));
    372:   PetscCall(DMDestroy(&dmM));
    373:   PetscCall(DMRestoreLocalVector(dm, &locU));
    374:   PetscCall(DMRestoreLocalVector(dmGrad, &locG));
    375:   PetscCall(DMDestroy(&dmGrad));

    377:   // TODO remove?
    378:   PetscScalar scale = 10.0;
    379:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-metric_scale", &scale, NULL));
    380:   PetscCall(VecScale(metric, scale));
    381:   PetscCall(VecViewFromOptions(metric, NULL, "-metric_view"));

    383:   DMLabel label = NULL;
    384:   PetscCall(DMGetLabel(dm, "marker", &label));
    385:   PetscCall(DMAdaptMetric(dm, metric, label, NULL, newdm));
    386:   PetscCall(VecDestroy(&metric));

    388:   PetscCall(DMViewFromOptions(*newdm, NULL, "-remesh_dmout_view"));

    390:   AppCtx *ctx;
    391:   PetscCall(DMGetApplicationContext(dm, &ctx));
    392:   PetscCall(DMSetApplicationContext(*newdm, ctx));
    393:   PetscCall(SetupDiscretization(*newdm, ctx));

    395:   // TODO
    396:   ((DM_Plex *)(*newdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
    397:   PetscFunctionReturn(PETSC_SUCCESS);
    398: }

    400: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
    401: {
    402:   DM        dm;
    403:   PetscReal t;

    405:   PetscFunctionBeginUser;
    406:   PetscCall(TSGetDM(ts, &dm));
    407:   PetscCall(TSGetTime(ts, &t));
    408:   PetscCall(DMComputeExactSolution(dm, t, u, NULL));
    409:   PetscCall(VecSetOptionsPrefix(u, NULL));
    410:   PetscFunctionReturn(PETSC_SUCCESS);
    411: }

    413: static PetscErrorCode TransferSetUp(TS ts, PetscInt step, PetscReal time, Vec U, PetscBool *remesh, void *tctx)
    414: {
    415:   AppCtx *ctx = (AppCtx *)tctx;

    417:   PetscFunctionBeginUser;
    418:   *remesh = PETSC_FALSE;
    419:   if (ctx->remesh_every > 0 && step && step % ctx->remesh_every == 0) {
    420:     DM dm;

    422:     *remesh = PETSC_TRUE;
    423:     PetscCall(TSGetDM(ts, &dm));
    424:     PetscCall(Remesh(dm, U, &ctx->remesh_dm));
    425:   }
    426:   PetscFunctionReturn(PETSC_SUCCESS);
    427: }

    429: static PetscErrorCode TransferVecs(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *tctx)
    430: {
    431:   AppCtx *ctx = (AppCtx *)tctx;
    432:   DM      dm;
    433:   Mat     Interp;

    435:   PetscFunctionBeginUser;
    436:   PetscCall(TSGetDM(ts, &dm));
    437:   PetscCall(DMCreateInterpolation(dm, ctx->remesh_dm, &Interp, NULL));
    438:   for (PetscInt i = 0; i < nv; i++) {
    439:     PetscCall(DMCreateGlobalVector(ctx->remesh_dm, &vout[i]));
    440:     PetscCall(MatMult(Interp, vin[i], vout[i]));
    441:   }
    442:   PetscCall(MatDestroy(&Interp));
    443:   PetscCall(TSSetDM(ts, ctx->remesh_dm));
    444:   PetscCall(DMDestroy(&ctx->remesh_dm));
    445:   PetscFunctionReturn(PETSC_SUCCESS);
    446: }

    448: int main(int argc, char **argv)
    449: {
    450:   DM     dm;
    451:   TS     ts;
    452:   Vec    u;
    453:   AppCtx ctx;

    455:   PetscFunctionBeginUser;
    456:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    457:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
    458:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm, &ctx));
    459:   PetscCall(DMSetApplicationContext(dm, &ctx));
    460:   PetscCall(SetupDiscretization(dm, &ctx));

    462:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    463:   PetscCall(TSSetDM(ts, dm));
    464:   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
    465:   if (ctx.expTS) {
    466:     PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFEM, &ctx));
    467:     if (ctx.lumped) PetscCall(DMTSCreateRHSMassMatrixLumped(dm));
    468:     else PetscCall(DMTSCreateRHSMassMatrix(dm));
    469:   } else {
    470:     PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
    471:     PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
    472:   }
    473:   PetscCall(TSSetMaxTime(ts, 1.0));
    474:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    475:   PetscCall(TSSetFromOptions(ts));
    476:   PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions));
    477:   PetscCall(TSSetResize(ts, PETSC_FALSE, TransferSetUp, TransferVecs, &ctx));

    479:   PetscCall(DMCreateGlobalVector(dm, &u));
    480:   PetscCall(DMTSCheckFromOptions(ts, u));
    481:   PetscCall(SetInitialConditions(ts, u));
    482:   PetscCall(PetscObjectSetName((PetscObject)u, "temperature"));

    484:   PetscCall(TSSetSolution(ts, u));
    485:   PetscCall(VecViewFromOptions(u, NULL, "-u0_view"));
    486:   PetscCall(VecDestroy(&u));
    487:   PetscCall(TSSolve(ts, NULL));

    489:   PetscCall(TSGetSolution(ts, &u));
    490:   PetscCall(VecViewFromOptions(u, NULL, "-uf_view"));
    491:   PetscCall(DMTSCheckFromOptions(ts, u));
    492:   if (ctx.expTS) PetscCall(DMTSDestroyRHSMassMatrix(dm));

    494:   PetscCall(TSDestroy(&ts));
    495:   PetscCall(DMDestroy(&dm));
    496:   PetscCall(PetscFinalize());
    497:   return 0;
    498: }

    500: /*TEST

    502:   test:
    503:     suffix: 2d_p1
    504:     requires: triangle
    505:     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
    506:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    507:   test:
    508:     suffix: 2d_p1_exp
    509:     requires: triangle
    510:     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -explicit \
    511:           -ts_type euler -ts_max_steps 4 -ts_time_step 1e-3 -ts_monitor_error
    512:   test:
    513:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
    514:     suffix: 2d_p1_sconv
    515:     requires: triangle
    516:     args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    517:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00001 -snes_error_if_not_converged -pc_type lu
    518:   test:
    519:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.1]
    520:     suffix: 2d_p1_sconv_2
    521:     requires: triangle
    522:     args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    523:           -ts_type beuler -ts_max_steps 1 -ts_time_step 1e-6 -snes_error_if_not_converged -pc_type lu
    524:   test:
    525:     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
    526:     suffix: 2d_p1_tconv
    527:     requires: triangle
    528:     args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
    529:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    530:   test:
    531:     # -dm_refine 6 -convest_num_refine 3 get L_2 convergence rate: [1.0]
    532:     suffix: 2d_p1_tconv_2
    533:     requires: triangle
    534:     args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
    535:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    536:   test:
    537:     # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid
    538:     suffix: 2d_p1_exp_tconv_2
    539:     requires: triangle
    540:     args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \
    541:           -ts_type euler -ts_max_steps 4 -ts_time_step 1e-4 -lumped 0 -mass_pc_type lu
    542:   test:
    543:     # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid
    544:     suffix: 2d_p1_exp_tconv_2_lump
    545:     requires: triangle
    546:     args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \
    547:           -ts_type euler -ts_max_steps 4 -ts_time_step 1e-4
    548:   test:
    549:     suffix: 2d_p2
    550:     requires: triangle
    551:     args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
    552:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    553:   test:
    554:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
    555:     suffix: 2d_p2_sconv
    556:     requires: triangle
    557:     args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    558:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00000001 -snes_error_if_not_converged -pc_type lu
    559:   test:
    560:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [3.1]
    561:     suffix: 2d_p2_sconv_2
    562:     requires: triangle
    563:     args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    564:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00000001 -snes_error_if_not_converged -pc_type lu
    565:   test:
    566:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
    567:     suffix: 2d_p2_tconv
    568:     requires: triangle
    569:     args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
    570:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    571:   test:
    572:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
    573:     suffix: 2d_p2_tconv_2
    574:     requires: triangle
    575:     args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
    576:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    577:   test:
    578:     suffix: 2d_q1
    579:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
    580:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    581:   test:
    582:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
    583:     suffix: 2d_q1_sconv
    584:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    585:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00001 -snes_error_if_not_converged -pc_type lu
    586:   test:
    587:     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
    588:     suffix: 2d_q1_tconv
    589:     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
    590:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    591:   test:
    592:     suffix: 2d_q2
    593:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
    594:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    595:   test:
    596:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
    597:     suffix: 2d_q2_sconv
    598:     args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    599:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00000001 -snes_error_if_not_converged -pc_type lu
    600:   test:
    601:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
    602:     suffix: 2d_q2_tconv
    603:     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
    604:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu

    606:   test:
    607:     suffix: 3d_p1
    608:     requires: ctetgen
    609:     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
    610:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    611:   test:
    612:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
    613:     suffix: 3d_p1_sconv
    614:     requires: ctetgen
    615:     args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    616:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00001 -snes_error_if_not_converged -pc_type lu
    617:   test:
    618:     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
    619:     suffix: 3d_p1_tconv
    620:     requires: ctetgen
    621:     args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
    622:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    623:   test:
    624:     suffix: 3d_p2
    625:     requires: ctetgen
    626:     args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
    627:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    628:   test:
    629:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
    630:     suffix: 3d_p2_sconv
    631:     requires: ctetgen
    632:     args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    633:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00000001 -snes_error_if_not_converged -pc_type lu
    634:   test:
    635:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
    636:     suffix: 3d_p2_tconv
    637:     requires: ctetgen
    638:     args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
    639:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    640:   test:
    641:     suffix: 3d_q1
    642:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
    643:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    644:   test:
    645:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
    646:     suffix: 3d_q1_sconv
    647:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    648:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00001 -snes_error_if_not_converged -pc_type lu
    649:   test:
    650:     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
    651:     suffix: 3d_q1_tconv
    652:     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
    653:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    654:   test:
    655:     suffix: 3d_q2
    656:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
    657:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu
    658:   test:
    659:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
    660:     suffix: 3d_q2_sconv
    661:     args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    662:           -ts_type beuler -ts_max_steps 1 -ts_time_step 0.00000001 -snes_error_if_not_converged -pc_type lu
    663:   test:
    664:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
    665:     suffix: 3d_q2_tconv
    666:     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
    667:           -ts_type beuler -ts_max_steps 4 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu

    669:   test:
    670:     # For a nice picture, -bd_dm_refine 2 -dm_refine 1 -dm_view hdf5:${PETSC_DIR}/sol.h5 -ts_monitor_solution hdf5:${PETSC_DIR}/sol.h5::append
    671:     suffix: egads_sphere
    672:     requires: egads ctetgen datafilespath
    673:     args: -sol_type quadratic_linear \
    674:           -dm_plex_boundary_filename ${DATAFILESPATH}/meshes/cad/sphere_example.egadslite -dm_plex_boundary_label marker \
    675:           -temp_petscspace_degree 2 -dmts_check .0001 \
    676:           -ts_type beuler -ts_max_steps 5 -ts_time_step 0.1 -snes_error_if_not_converged -pc_type lu

    678:   test:
    679:     suffix: remesh
    680:     requires: triangle mmg
    681:     args: -sol_type quadratic_trig -dm_refine 2 -temp_petscspace_degree 1 -ts_type beuler -ts_time_step 0.01 -snes_error_if_not_converged -pc_type lu -grad_petscspace_degree 1 -dm_adaptor mmg -dm_plex_hash_location -remesh_every 5
    682:     output_file: output/empty.out

    684: TEST*/