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