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 {SOL_QUADRATIC_LINEAR, SOL_QUADRATIC_TRIG, SOL_TRIG_LINEAR, SOL_TRIG_TRIG, NUM_SOLUTION_TYPES} SolutionType;
 17: const char *solutionTypes[NUM_SOLUTION_TYPES+1] = {"quadratic_linear", "quadratic_trig", "trig_linear", "trig_trig", "unknown"};

 19: typedef struct {
 20:   SolutionType solType; /* Type of exact solution */
 21:   /* Solver setup */
 22:   PetscBool expTS;      /* Flag for explicit timestepping */
 23:   PetscBool lumped;     /* Lump the mass matrix */
 24: } AppCtx;

 26: /*
 27: Exact 2D solution:
 28:   u    = 2t + x^2 + y^2
 29:   u_t  = 2
 30:   \Delta u = 2 + 2 = 4
 31:   f    = 2
 32:   F(u) = 2 - (2 + 2) + 2 = 0

 34: Exact 3D solution:
 35:   u = 3t + x^2 + y^2 + z^2
 36:   F(u) = 3 - (2 + 2 + 2) + 3 = 0
 37: */
 38: static PetscErrorCode mms_quad_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 39: {
 40:   PetscInt d;

 42:   *u = dim*time;
 43:   for (d = 0; d < dim; ++d) *u += x[d]*x[d];
 44:   return 0;
 45: }

 47: static PetscErrorCode mms_quad_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 48: {
 49:   *u = dim;
 50:   return 0;
 51: }

 53: static void f0_quad_lin_exp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 54:                             const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 55:                             const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 56:                             PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 57: {
 58:   f0[0] = -(PetscScalar) dim;
 59: }
 60: static void f0_quad_lin(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 61:                         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 62:                         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 63:                         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 64: {
 65:   PetscScalar exp[1] = {0.};
 66:   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);
 67:   f0[0] = u_t[0] - exp[0];
 68: }

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

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

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

 88: static PetscErrorCode mms_quad_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 89: {
 90:   *u = -dim*PetscSinReal(time);
 91:   return 0;
 92: }

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

111: /*
112: Exact 2D solution:
113:   u = 2\pi^2 t + cos(\pi x) + cos(\pi y)
114:   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

116: Exact 3D solution:
117:   u = 3\pi^2 t + cos(\pi x) + cos(\pi y) + cos(\pi z)
118:   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
119: */
120: static PetscErrorCode mms_trig_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
121: {
122:   PetscInt d;

124:   *u = dim*PetscSqr(PETSC_PI)*time;
125:   for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI*x[d]);
126:   return 0;
127: }

129: static PetscErrorCode mms_trig_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
130: {
131:   *u = dim*PetscSqr(PETSC_PI);
132:   return 0;
133: }

135: static void f0_trig_lin(PetscInt dim, PetscInt Nf, PetscInt NfAux,
136:                         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
137:                         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
138:                         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
139: {
140:   PetscInt d;
141:   f0[0] = u_t[0];
142:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI)*(PetscCosReal(PETSC_PI*x[d]) - 1.0);
143: }

145: /*
146: Exact 2D solution:
147:   u    = pi^2 cos(t) + cos(\pi x) + cos(\pi y)
148:   u_t  = -pi^2 sin(t)
149:   \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y))
150:   f    = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y))
151:   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

153: Exact 3D solution:
154:   u    = pi^2 cos(t) + cos(\pi x) + cos(\pi y) + cos(\pi z)
155:   u_t  = -pi^2 sin(t)
156:   \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z))
157:   f    = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z))
158:   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
159: */
160: static PetscErrorCode mms_trig_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
161: {
162:   PetscInt d;

164:   *u = PetscSqr(PETSC_PI)*PetscCosReal(time);
165:   for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI*x[d]);
166:   return 0;
167: }

169: static PetscErrorCode mms_trig_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
170: {
171:   *u = -PetscSqr(PETSC_PI)*PetscSinReal(time);
172:   return 0;
173: }

175: static void f0_trig_trig_exp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
176:                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
177:                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
178:                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
179: {
180:   PetscInt d;
181:   f0[0] -= PetscSqr(PETSC_PI)*PetscSinReal(t);
182:   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI)*PetscCosReal(PETSC_PI*x[d]);
183: }
184: static void f0_trig_trig(PetscInt dim, PetscInt Nf, PetscInt NfAux,
185:                          const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
186:                          const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
187:                          PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
188: {
189:   PetscScalar exp[1] = {0.};
190:   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);
191:   f0[0] = u_t[0] - exp[0];
192: }

194: static void f1_temp_exp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
195:                         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
196:                         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
197:                         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
198: {
199:   PetscInt d;
200:   for (d = 0; d < dim; ++d) f1[d] = -u_x[d];
201: }
202: static void f1_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
203:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
204:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
205:                     PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
206: {
207:   PetscInt d;
208:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
209: }

211: static void g3_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
212:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
213:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
214:                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
215: {
216:   PetscInt d;
217:   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
218: }

220: static void g0_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
221:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
222:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
223:                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
224: {
225:   g0[0] = u_tShift*1.0;
226: }

228: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
229: {
230:   PetscInt       sol;

234:   options->solType  = SOL_QUADRATIC_LINEAR;
235:   options->expTS    = PETSC_FALSE;
236:   options->lumped   = PETSC_TRUE;

238:   PetscOptionsBegin(comm, "", "Heat Equation Options", "DMPLEX");
239:   PetscOptionsEList("-sol_type", "Type of exact solution", "ex45.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);
240:   options->solType = (SolutionType) sol;
241:   PetscOptionsBool("-explicit", "Use explicit timestepping", "ex45.c", options->expTS, &options->expTS, NULL);
242:   PetscOptionsBool("-lumped", "Lump the mass matrix", "ex45.c", options->lumped, &options->lumped, NULL);
243:   PetscOptionsEnd();
244:   return 0;
245: }

247: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
248: {
250:   DMCreate(comm, dm);
251:   DMSetType(*dm, DMPLEX);
252:   DMSetFromOptions(*dm);
253:   DMViewFromOptions(*dm, NULL, "-dm_view");
254:   return 0;
255: }

257: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
258: {
259:   PetscDS        ds;
260:   DMLabel        label;
261:   const PetscInt id = 1;

264:   DMGetLabel(dm, "marker", &label);
265:   DMGetDS(dm, &ds);
266:   PetscDSSetJacobian(ds, 0, 0, g0_temp, NULL, NULL, g3_temp);
267:   switch (ctx->solType) {
268:     case SOL_QUADRATIC_LINEAR:
269:       PetscDSSetResidual(ds, 0, f0_quad_lin,  f1_temp);
270:       PetscDSSetRHSResidual(ds, 0, f0_quad_lin_exp, f1_temp_exp);
271:       PetscDSSetExactSolution(ds, 0, mms_quad_lin, ctx);
272:       PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_lin_t, ctx);
273:       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);
274:       break;
275:     case SOL_QUADRATIC_TRIG:
276:       PetscDSSetResidual(ds, 0, f0_quad_trig, f1_temp);
277:       PetscDSSetRHSResidual(ds, 0, f0_quad_trig_exp, f1_temp_exp);
278:       PetscDSSetExactSolution(ds, 0, mms_quad_trig, ctx);
279:       PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_trig_t, ctx);
280:       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);
281:       break;
282:     case SOL_TRIG_LINEAR:
283:       PetscDSSetResidual(ds, 0, f0_trig_lin,  f1_temp);
284:       PetscDSSetExactSolution(ds, 0, mms_trig_lin, ctx);
285:       PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_lin_t, ctx);
286:       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);
287:       break;
288:     case SOL_TRIG_TRIG:
289:       PetscDSSetResidual(ds, 0, f0_trig_trig, f1_temp);
290:       PetscDSSetRHSResidual(ds, 0, f0_trig_trig_exp, f1_temp_exp);
291:       PetscDSSetExactSolution(ds, 0, mms_trig_trig, ctx);
292:       PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_trig_t, ctx);
293:       break;
294:     default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType);
295:   }
296:   return 0;
297: }

299: static PetscErrorCode SetupDiscretization(DM dm, AppCtx* ctx)
300: {
301:   DM             cdm = dm;
302:   PetscFE        fe;
303:   DMPolytopeType ct;
304:   PetscBool      simplex;
305:   PetscInt       dim, cStart;

308:   DMGetDimension(dm, &dim);
309:   DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
310:   DMPlexGetCellType(dm, cStart, &ct);
311:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
312:   /* Create finite element */
313:   PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", -1, &fe);
314:   PetscObjectSetName((PetscObject) fe, "temperature");
315:   /* Set discretization and boundary conditions for each mesh */
316:   DMSetField(dm, 0, NULL, (PetscObject) fe);
317:   DMCreateDS(dm);
318:   if (ctx->expTS) {
319:     PetscDS ds;

321:     DMGetDS(dm, &ds);
322:     PetscDSSetImplicit(ds, 0, PETSC_FALSE);
323:   }
324:   SetupProblem(dm, ctx);
325:   while (cdm) {
326:     DMCopyDisc(dm, cdm);
327:     DMGetCoarseDM(cdm, &cdm);
328:   }
329:   PetscFEDestroy(&fe);
330:   return 0;
331: }

333: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
334: {
335:   DM             dm;
336:   PetscReal      t;

338:   TSGetDM(ts, &dm);
339:   TSGetTime(ts, &t);
340:   DMComputeExactSolution(dm, t, u, NULL);
341:   return 0;
342: }

344: int main(int argc, char **argv)
345: {
346:   DM             dm;
347:   TS             ts;
348:   Vec            u;
349:   AppCtx         ctx;

351:   PetscInitialize(&argc, &argv, NULL, help);
352:   ProcessOptions(PETSC_COMM_WORLD, &ctx);
353:   CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);
354:   DMSetApplicationContext(dm, &ctx);
355:   SetupDiscretization(dm, &ctx);

357:   TSCreate(PETSC_COMM_WORLD, &ts);
358:   TSSetDM(ts, dm);
359:   DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
360:   if (ctx.expTS) {
361:     DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFEM, &ctx);
362:     if (ctx.lumped) DMTSCreateRHSMassMatrixLumped(dm);
363:     else            DMTSCreateRHSMassMatrix(dm);
364:   } else {
365:     DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
366:     DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
367:   }
368:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
369:   TSSetFromOptions(ts);
370:   TSSetComputeInitialCondition(ts, SetInitialConditions);

372:   DMCreateGlobalVector(dm, &u);
373:   DMTSCheckFromOptions(ts, u);
374:   SetInitialConditions(ts, u);
375:   PetscObjectSetName((PetscObject) u, "temperature");
376:   TSSolve(ts, u);
377:   DMTSCheckFromOptions(ts, u);
378:   if (ctx.expTS) DMTSDestroyRHSMassMatrix(dm);

380:   VecDestroy(&u);
381:   TSDestroy(&ts);
382:   DMDestroy(&dm);
383:   PetscFinalize();
384:   return 0;
385: }

387: /*TEST

389:   test:
390:     suffix: 2d_p1
391:     requires: triangle
392:     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
393:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
394:   test:
395:     suffix: 2d_p1_exp
396:     requires: triangle
397:     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -explicit \
398:           -ts_type euler -ts_max_steps 4 -ts_dt 1e-3 -ts_monitor_error
399:   test:
400:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
401:     suffix: 2d_p1_sconv
402:     requires: triangle
403:     args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
404:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
405:   test:
406:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.1]
407:     suffix: 2d_p1_sconv_2
408:     requires: triangle
409:     args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
410:           -ts_type beuler -ts_max_steps 1 -ts_dt 1e-6 -snes_error_if_not_converged -pc_type lu
411:   test:
412:     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
413:     suffix: 2d_p1_tconv
414:     requires: triangle
415:     args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
416:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
417:   test:
418:     # -dm_refine 6 -convest_num_refine 3 get L_2 convergence rate: [1.0]
419:     suffix: 2d_p1_tconv_2
420:     requires: triangle
421:     args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
422:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
423:   test:
424:     # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid
425:     suffix: 2d_p1_exp_tconv_2
426:     requires: triangle
427:     args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \
428:           -ts_type euler -ts_max_steps 4 -ts_dt 1e-4 -lumped 0 -mass_pc_type lu
429:   test:
430:     # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid
431:     suffix: 2d_p1_exp_tconv_2_lump
432:     requires: triangle
433:     args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \
434:           -ts_type euler -ts_max_steps 4 -ts_dt 1e-4
435:   test:
436:     suffix: 2d_p2
437:     requires: triangle
438:     args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
439:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
440:   test:
441:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
442:     suffix: 2d_p2_sconv
443:     requires: triangle
444:     args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
445:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
446:   test:
447:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [3.1]
448:     suffix: 2d_p2_sconv_2
449:     requires: triangle
450:     args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
451:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
452:   test:
453:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
454:     suffix: 2d_p2_tconv
455:     requires: triangle
456:     args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
457:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
458:   test:
459:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
460:     suffix: 2d_p2_tconv_2
461:     requires: triangle
462:     args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
463:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
464:   test:
465:     suffix: 2d_q1
466:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
467:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
468:   test:
469:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
470:     suffix: 2d_q1_sconv
471:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
472:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
473:   test:
474:     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
475:     suffix: 2d_q1_tconv
476:     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
477:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
478:   test:
479:     suffix: 2d_q2
480:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
481:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
482:   test:
483:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
484:     suffix: 2d_q2_sconv
485:     args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
486:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
487:   test:
488:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
489:     suffix: 2d_q2_tconv
490:     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
491:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu

493:   test:
494:     suffix: 3d_p1
495:     requires: ctetgen
496:     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
497:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
498:   test:
499:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
500:     suffix: 3d_p1_sconv
501:     requires: ctetgen
502:     args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
503:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
504:   test:
505:     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
506:     suffix: 3d_p1_tconv
507:     requires: ctetgen
508:     args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
509:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
510:   test:
511:     suffix: 3d_p2
512:     requires: ctetgen
513:     args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
514:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
515:   test:
516:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
517:     suffix: 3d_p2_sconv
518:     requires: ctetgen
519:     args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
520:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
521:   test:
522:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
523:     suffix: 3d_p2_tconv
524:     requires: ctetgen
525:     args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
526:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
527:   test:
528:     suffix: 3d_q1
529:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
530:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
531:   test:
532:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
533:     suffix: 3d_q1_sconv
534:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
535:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
536:   test:
537:     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
538:     suffix: 3d_q1_tconv
539:     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
540:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
541:   test:
542:     suffix: 3d_q2
543:     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
544:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
545:   test:
546:     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
547:     suffix: 3d_q2_sconv
548:     args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
549:           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
550:   test:
551:     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
552:     suffix: 3d_q2_tconv
553:     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
554:           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu

556:   test:
557:     # 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
558:     suffix: egads_sphere
559:     requires: egads ctetgen
560:     args: -sol_type quadratic_linear \
561:           -dm_plex_boundary_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_plex_boundary_label marker -bd_dm_plex_scale 40 \
562:           -temp_petscspace_degree 2 -dmts_check .0001 \
563:           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu

565: TEST*/