xref: /petsc/src/ts/tutorials/ex45.c (revision 9371c9d470a9602b6d10a8bf50c9b2280a79e45a)
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";
5 
6 #include <petscdmplex.h>
7 #include <petscds.h>
8 #include <petscts.h>
9 
10 /*
11   Heat equation:
12 
13     du/dt - \Delta u + f = 0
14 */
15 
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"};
24 
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 } AppCtx;
31 
32 /*
33 Exact 2D solution:
34   u    = 2t + x^2 + y^2
35   u_t  = 2
36   \Delta u = 2 + 2 = 4
37   f    = 2
38   F(u) = 2 - (2 + 2) + 2 = 0
39 
40 Exact 3D solution:
41   u = 3t + x^2 + y^2 + z^2
42   F(u) = 3 - (2 + 2 + 2) + 3 = 0
43 */
44 static PetscErrorCode mms_quad_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) {
45   PetscInt d;
46 
47   *u = dim * time;
48   for (d = 0; d < dim; ++d) *u += x[d] * x[d];
49   return 0;
50 }
51 
52 static PetscErrorCode mms_quad_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) {
53   *u = dim;
54   return 0;
55 }
56 
57 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[]) {
58   f0[0] = -(PetscScalar)dim;
59 }
60 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[]) {
61   PetscScalar exp[1] = {0.};
62   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);
63   f0[0] = u_t[0] - exp[0];
64 }
65 
66 /*
67 Exact 2D solution:
68   u = 2*cos(t) + x^2 + y^2
69   F(u) = -2*sint(t) - (2 + 2) + 2*sin(t) + 4 = 0
70 
71 Exact 3D solution:
72   u = 3*cos(t) + x^2 + y^2 + z^2
73   F(u) = -3*sin(t) - (2 + 2 + 2) + 3*sin(t) + 6 = 0
74 */
75 static PetscErrorCode mms_quad_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) {
76   PetscInt d;
77 
78   *u = dim * PetscCosReal(time);
79   for (d = 0; d < dim; ++d) *u += x[d] * x[d];
80   return 0;
81 }
82 
83 static PetscErrorCode mms_quad_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) {
84   *u = -dim * PetscSinReal(time);
85   return 0;
86 }
87 
88 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[]) {
89   f0[0] = -dim * (PetscSinReal(t) + 2.0);
90 }
91 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[]) {
92   PetscScalar exp[1] = {0.};
93   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);
94   f0[0] = u_t[0] - exp[0];
95 }
96 
97 /*
98 Exact 2D solution:
99   u = 2\pi^2 t + cos(\pi x) + cos(\pi y)
100   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
101 
102 Exact 3D solution:
103   u = 3\pi^2 t + cos(\pi x) + cos(\pi y) + cos(\pi z)
104   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
105 */
106 static PetscErrorCode mms_trig_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) {
107   PetscInt d;
108 
109   *u = dim * PetscSqr(PETSC_PI) * time;
110   for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI * x[d]);
111   return 0;
112 }
113 
114 static PetscErrorCode mms_trig_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) {
115   *u = dim * PetscSqr(PETSC_PI);
116   return 0;
117 }
118 
119 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[]) {
120   PetscInt d;
121   f0[0] = u_t[0];
122   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI) * (PetscCosReal(PETSC_PI * x[d]) - 1.0);
123 }
124 
125 /*
126 Exact 2D solution:
127   u    = pi^2 cos(t) + cos(\pi x) + cos(\pi y)
128   u_t  = -pi^2 sin(t)
129   \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y))
130   f    = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y))
131   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
132 
133 Exact 3D solution:
134   u    = pi^2 cos(t) + cos(\pi x) + cos(\pi y) + cos(\pi z)
135   u_t  = -pi^2 sin(t)
136   \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z))
137   f    = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z))
138   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
139 */
140 static PetscErrorCode mms_trig_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) {
141   PetscInt d;
142 
143   *u = PetscSqr(PETSC_PI) * PetscCosReal(time);
144   for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI * x[d]);
145   return 0;
146 }
147 
148 static PetscErrorCode mms_trig_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) {
149   *u = -PetscSqr(PETSC_PI) * PetscSinReal(time);
150   return 0;
151 }
152 
153 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[]) {
154   PetscInt d;
155   f0[0] -= PetscSqr(PETSC_PI) * PetscSinReal(t);
156   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI) * PetscCosReal(PETSC_PI * x[d]);
157 }
158 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[]) {
159   PetscScalar exp[1] = {0.};
160   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);
161   f0[0] = u_t[0] - exp[0];
162 }
163 
164 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[]) {
165   PetscInt d;
166   for (d = 0; d < dim; ++d) f1[d] = -u_x[d];
167 }
168 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[]) {
169   PetscInt d;
170   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
171 }
172 
173 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[]) {
174   PetscInt d;
175   for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
176 }
177 
178 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[]) {
179   g0[0] = u_tShift * 1.0;
180 }
181 
182 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) {
183   PetscInt sol;
184 
185   PetscFunctionBeginUser;
186   options->solType = SOL_QUADRATIC_LINEAR;
187   options->expTS   = PETSC_FALSE;
188   options->lumped  = PETSC_TRUE;
189 
190   PetscOptionsBegin(comm, "", "Heat Equation Options", "DMPLEX");
191   PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex45.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL));
192   options->solType = (SolutionType)sol;
193   PetscCall(PetscOptionsBool("-explicit", "Use explicit timestepping", "ex45.c", options->expTS, &options->expTS, NULL));
194   PetscCall(PetscOptionsBool("-lumped", "Lump the mass matrix", "ex45.c", options->lumped, &options->lumped, NULL));
195   PetscOptionsEnd();
196   PetscFunctionReturn(0);
197 }
198 
199 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx) {
200   PetscFunctionBeginUser;
201   PetscCall(DMCreate(comm, dm));
202   PetscCall(DMSetType(*dm, DMPLEX));
203   PetscCall(DMSetFromOptions(*dm));
204   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
205   PetscFunctionReturn(0);
206 }
207 
208 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) {
209   PetscDS        ds;
210   DMLabel        label;
211   const PetscInt id = 1;
212 
213   PetscFunctionBeginUser;
214   PetscCall(DMGetLabel(dm, "marker", &label));
215   PetscCall(DMGetDS(dm, &ds));
216   PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_temp, NULL, NULL, g3_temp));
217   switch (ctx->solType) {
218   case SOL_QUADRATIC_LINEAR:
219     PetscCall(PetscDSSetResidual(ds, 0, f0_quad_lin, f1_temp));
220     PetscCall(PetscDSSetRHSResidual(ds, 0, f0_quad_lin_exp, f1_temp_exp));
221     PetscCall(PetscDSSetExactSolution(ds, 0, mms_quad_lin, ctx));
222     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_lin_t, ctx));
223     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));
224     break;
225   case SOL_QUADRATIC_TRIG:
226     PetscCall(PetscDSSetResidual(ds, 0, f0_quad_trig, f1_temp));
227     PetscCall(PetscDSSetRHSResidual(ds, 0, f0_quad_trig_exp, f1_temp_exp));
228     PetscCall(PetscDSSetExactSolution(ds, 0, mms_quad_trig, ctx));
229     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_trig_t, ctx));
230     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));
231     break;
232   case SOL_TRIG_LINEAR:
233     PetscCall(PetscDSSetResidual(ds, 0, f0_trig_lin, f1_temp));
234     PetscCall(PetscDSSetExactSolution(ds, 0, mms_trig_lin, ctx));
235     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_lin_t, ctx));
236     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));
237     break;
238   case SOL_TRIG_TRIG:
239     PetscCall(PetscDSSetResidual(ds, 0, f0_trig_trig, f1_temp));
240     PetscCall(PetscDSSetRHSResidual(ds, 0, f0_trig_trig_exp, f1_temp_exp));
241     PetscCall(PetscDSSetExactSolution(ds, 0, mms_trig_trig, ctx));
242     PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_trig_t, ctx));
243     break;
244   default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType);
245   }
246   PetscFunctionReturn(0);
247 }
248 
249 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx) {
250   DM             cdm = dm;
251   PetscFE        fe;
252   DMPolytopeType ct;
253   PetscBool      simplex;
254   PetscInt       dim, cStart;
255 
256   PetscFunctionBeginUser;
257   PetscCall(DMGetDimension(dm, &dim));
258   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
259   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
260   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
261   /* Create finite element */
262   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", -1, &fe));
263   PetscCall(PetscObjectSetName((PetscObject)fe, "temperature"));
264   /* Set discretization and boundary conditions for each mesh */
265   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
266   PetscCall(DMCreateDS(dm));
267   if (ctx->expTS) {
268     PetscDS ds;
269 
270     PetscCall(DMGetDS(dm, &ds));
271     PetscCall(PetscDSSetImplicit(ds, 0, PETSC_FALSE));
272   }
273   PetscCall(SetupProblem(dm, ctx));
274   while (cdm) {
275     PetscCall(DMCopyDisc(dm, cdm));
276     PetscCall(DMGetCoarseDM(cdm, &cdm));
277   }
278   PetscCall(PetscFEDestroy(&fe));
279   PetscFunctionReturn(0);
280 }
281 
282 static PetscErrorCode SetInitialConditions(TS ts, Vec u) {
283   DM        dm;
284   PetscReal t;
285 
286   PetscFunctionBeginUser;
287   PetscCall(TSGetDM(ts, &dm));
288   PetscCall(TSGetTime(ts, &t));
289   PetscCall(DMComputeExactSolution(dm, t, u, NULL));
290   PetscFunctionReturn(0);
291 }
292 
293 int main(int argc, char **argv) {
294   DM     dm;
295   TS     ts;
296   Vec    u;
297   AppCtx ctx;
298 
299   PetscFunctionBeginUser;
300   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
301   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
302   PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm, &ctx));
303   PetscCall(DMSetApplicationContext(dm, &ctx));
304   PetscCall(SetupDiscretization(dm, &ctx));
305 
306   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
307   PetscCall(TSSetDM(ts, dm));
308   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
309   if (ctx.expTS) {
310     PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFEM, &ctx));
311     if (ctx.lumped) PetscCall(DMTSCreateRHSMassMatrixLumped(dm));
312     else PetscCall(DMTSCreateRHSMassMatrix(dm));
313   } else {
314     PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
315     PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
316   }
317   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
318   PetscCall(TSSetFromOptions(ts));
319   PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions));
320 
321   PetscCall(DMCreateGlobalVector(dm, &u));
322   PetscCall(DMTSCheckFromOptions(ts, u));
323   PetscCall(SetInitialConditions(ts, u));
324   PetscCall(PetscObjectSetName((PetscObject)u, "temperature"));
325   PetscCall(TSSolve(ts, u));
326   PetscCall(DMTSCheckFromOptions(ts, u));
327   if (ctx.expTS) PetscCall(DMTSDestroyRHSMassMatrix(dm));
328 
329   PetscCall(VecDestroy(&u));
330   PetscCall(TSDestroy(&ts));
331   PetscCall(DMDestroy(&dm));
332   PetscCall(PetscFinalize());
333   return 0;
334 }
335 
336 /*TEST
337 
338   test:
339     suffix: 2d_p1
340     requires: triangle
341     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
342           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
343   test:
344     suffix: 2d_p1_exp
345     requires: triangle
346     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -explicit \
347           -ts_type euler -ts_max_steps 4 -ts_dt 1e-3 -ts_monitor_error
348   test:
349     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
350     suffix: 2d_p1_sconv
351     requires: triangle
352     args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
353           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
354   test:
355     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.1]
356     suffix: 2d_p1_sconv_2
357     requires: triangle
358     args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
359           -ts_type beuler -ts_max_steps 1 -ts_dt 1e-6 -snes_error_if_not_converged -pc_type lu
360   test:
361     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
362     suffix: 2d_p1_tconv
363     requires: triangle
364     args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
365           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
366   test:
367     # -dm_refine 6 -convest_num_refine 3 get L_2 convergence rate: [1.0]
368     suffix: 2d_p1_tconv_2
369     requires: triangle
370     args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
371           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
372   test:
373     # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid
374     suffix: 2d_p1_exp_tconv_2
375     requires: triangle
376     args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \
377           -ts_type euler -ts_max_steps 4 -ts_dt 1e-4 -lumped 0 -mass_pc_type lu
378   test:
379     # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid
380     suffix: 2d_p1_exp_tconv_2_lump
381     requires: triangle
382     args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \
383           -ts_type euler -ts_max_steps 4 -ts_dt 1e-4
384   test:
385     suffix: 2d_p2
386     requires: triangle
387     args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
388           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
389   test:
390     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
391     suffix: 2d_p2_sconv
392     requires: triangle
393     args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
394           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
395   test:
396     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [3.1]
397     suffix: 2d_p2_sconv_2
398     requires: triangle
399     args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
400           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
401   test:
402     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
403     suffix: 2d_p2_tconv
404     requires: triangle
405     args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
406           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
407   test:
408     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
409     suffix: 2d_p2_tconv_2
410     requires: triangle
411     args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
412           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
413   test:
414     suffix: 2d_q1
415     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
416           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
417   test:
418     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
419     suffix: 2d_q1_sconv
420     args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
421           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
422   test:
423     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
424     suffix: 2d_q1_tconv
425     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
426           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
427   test:
428     suffix: 2d_q2
429     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
430           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
431   test:
432     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
433     suffix: 2d_q2_sconv
434     args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
435           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
436   test:
437     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
438     suffix: 2d_q2_tconv
439     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
440           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
441 
442   test:
443     suffix: 3d_p1
444     requires: ctetgen
445     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
446           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
447   test:
448     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
449     suffix: 3d_p1_sconv
450     requires: ctetgen
451     args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
452           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
453   test:
454     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
455     suffix: 3d_p1_tconv
456     requires: ctetgen
457     args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
458           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
459   test:
460     suffix: 3d_p2
461     requires: ctetgen
462     args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
463           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
464   test:
465     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
466     suffix: 3d_p2_sconv
467     requires: ctetgen
468     args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
469           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
470   test:
471     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
472     suffix: 3d_p2_tconv
473     requires: ctetgen
474     args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
475           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
476   test:
477     suffix: 3d_q1
478     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
479           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
480   test:
481     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
482     suffix: 3d_q1_sconv
483     args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
484           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
485   test:
486     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
487     suffix: 3d_q1_tconv
488     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
489           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
490   test:
491     suffix: 3d_q2
492     args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
493           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
494   test:
495     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
496     suffix: 3d_q2_sconv
497     args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
498           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
499   test:
500     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
501     suffix: 3d_q2_tconv
502     args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
503           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
504 
505   test:
506     # 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
507     suffix: egads_sphere
508     requires: egads ctetgen
509     args: -sol_type quadratic_linear \
510           -dm_plex_boundary_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_plex_boundary_label marker -bd_dm_plex_scale 40 \
511           -temp_petscspace_degree 2 -dmts_check .0001 \
512           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
513 
514 TEST*/
515