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, void *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, void *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, void *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, void *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, void *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, void *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, void *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, void *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, (void (*)(void))mms_quad_lin, (void (*)(void))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, (void (*)(void))mms_quad_trig, (void (*)(void))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, (void (*)(void))mms_trig_lin, (void (*)(void))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, 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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_dt 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
673:     args: -sol_type quadratic_linear \
674:           -dm_plex_boundary_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/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_dt 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_dt 0.01 -snes_error_if_not_converged -pc_type lu -grad_petscspace_degree 1 -dm_adaptor mmg -dm_plex_hash_location -remesh_every 5

683: TEST*/