xref: /petsc/src/ts/tutorials/ex45.c (revision 0bb0ac0dea362ee349c21557f7fd33dd750a14b5)
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 {SOL_QUADRATIC_LINEAR, SOL_QUADRATIC_TRIG, SOL_TRIG_LINEAR, NUM_SOLUTION_TYPES} SolutionType;
17 const char *solutionTypes[NUM_SOLUTION_TYPES+1] = {"quadratic_linear", "quadratic_trig", "trig_linear", "unknown"};
18 
19 typedef struct {
20   char         filename[PETSC_MAX_PATH_LEN];   /* Mesh filename */
21   char         bdfilename[PETSC_MAX_PATH_LEN]; /* Mesh boundary filename */
22   PetscReal    scale;                          /* Scale factor for mesh */
23   SolutionType solType;                        /* Type of exact solution */
24 } AppCtx;
25 
26 /*
27 Exact 2D solution:
28   u = 2t + x^2 + y^2
29   F(u) = 2 - (2 + 2) + 2 = 0
30 
31 Exact 3D solution:
32   u = 3t + x^2 + y^2 + z^2
33   F(u) = 3 - (2 + 2 + 2) + 3 = 0
34 */
35 static PetscErrorCode mms_quad_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
36 {
37   PetscInt d;
38 
39   *u = dim*time;
40   for (d = 0; d < dim; ++d) *u += x[d]*x[d];
41   return 0;
42 }
43 
44 static PetscErrorCode mms_quad_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
45 {
46   *u = dim;
47   return 0;
48 }
49 
50 static void f0_quad_lin(PetscInt dim, PetscInt Nf, PetscInt NfAux,
51                         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
52                         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
53                         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
54 {
55   f0[0] = u_t[0] + (PetscScalar) dim;
56 }
57 
58 /*
59 Exact 2D solution:
60   u = 2*cos(t) + x^2 + y^2
61   F(u) = -2*sint(t) - (2 + 2) + 2*sin(t) + 4 = 0
62 
63 Exact 3D solution:
64   u = 3*cos(t) + x^2 + y^2 + z^2
65   F(u) = -3*sin(t) - (2 + 2 + 2) + 3*sin(t) + 6 = 0
66 */
67 static PetscErrorCode mms_quad_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
68 {
69   PetscInt d;
70 
71   *u = dim*PetscCosReal(time);
72   for (d = 0; d < dim; ++d) *u += x[d]*x[d];
73   return 0;
74 }
75 
76 static PetscErrorCode mms_quad_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
77 {
78   *u = -dim*PetscSinReal(time);
79   return 0;
80 }
81 
82 static void f0_quad_trig(PetscInt dim, PetscInt Nf, PetscInt NfAux,
83                          const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
84                          const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
85                          PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
86 {
87   f0[0] = u_t[0] + dim*(PetscSinReal(t) + 2.0);
88 }
89 
90 /*
91 Exact 2D solution:
92   u = 2\pi^2 t + cos(\pi x) + cos(\pi y)
93   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
94 
95 Exact 3D solution:
96   u = 3\pi^2 t + cos(\pi x) + cos(\pi y) + cos(\pi z)
97   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
98 */
99 static PetscErrorCode mms_trig_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
100 {
101   PetscInt d;
102 
103   *u = dim*PetscSqr(PETSC_PI)*time;
104   for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI*x[d]);
105   return 0;
106 }
107 
108 static PetscErrorCode mms_trig_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
109 {
110   *u = dim*PetscSqr(PETSC_PI);
111   return 0;
112 }
113 
114 static void f0_trig_lin(PetscInt dim, PetscInt Nf, PetscInt NfAux,
115                         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
116                         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
117                         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
118 {
119   PetscInt d;
120   f0[0] = u_t[0];
121   for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI)*(PetscCosReal(PETSC_PI*x[d]) - 1.0);
122 }
123 
124 static void f1_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
125                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
126                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
127                     PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
128 {
129   PetscInt d;
130   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
131 }
132 
133 static void g3_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
134                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
135                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
136                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
137 {
138   PetscInt d;
139   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
140 }
141 
142 static void g0_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
143                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
144                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
145                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
146 {
147   g0[0] = u_tShift*1.0;
148 }
149 
150 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
151 {
152   PetscInt       sol;
153   PetscErrorCode ierr;
154 
155   PetscFunctionBeginUser;
156   options->filename[0]   = '\0';
157   options->bdfilename[0] = '\0';
158   options->scale         = 0.0;
159   options->solType       = SOL_QUADRATIC_LINEAR;
160 
161   ierr = PetscOptionsBegin(comm, "", "Heat Equation Options", "DMPLEX");CHKERRQ(ierr);
162   ierr = PetscOptionsString("-filename", "The mesh file", "ex45.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr);
163   ierr = PetscOptionsString("-bd_filename", "The mesh boundary file", "ex45.c", options->bdfilename, options->bdfilename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr);
164   ierr = PetscOptionsReal("-scale", "Scale factor for the mesh", "ex45.c", options->scale, &options->scale, NULL);CHKERRQ(ierr);
165   sol  = options->solType;
166   ierr = PetscOptionsEList("-sol_type", "Type of exact solution", "ex45.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);CHKERRQ(ierr);
167   options->solType = (SolutionType) sol;
168   ierr = PetscOptionsEnd();CHKERRQ(ierr);
169   PetscFunctionReturn(0);
170 }
171 
172 static PetscErrorCode CreateBCLabel(DM dm, const char name[])
173 {
174   DM             plex;
175   DMLabel        label;
176   PetscBool      hasLabel;
177   PetscErrorCode ierr;
178 
179   PetscFunctionBeginUser;
180   ierr = DMHasLabel(dm, name, &hasLabel);CHKERRQ(ierr);
181   if (hasLabel) PetscFunctionReturn(0);
182   ierr = DMCreateLabel(dm, name);CHKERRQ(ierr);
183   ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr);
184   ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr);
185   ierr = DMPlexMarkBoundaryFaces(plex, 1, label);CHKERRQ(ierr);
186   ierr = DMDestroy(&plex);CHKERRQ(ierr);
187   PetscFunctionReturn(0);
188 }
189 
190 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
191 {
192   size_t         len, lenbd;
193   PetscErrorCode ierr;
194 
195   PetscFunctionBeginUser;
196   ierr = PetscStrlen(ctx->filename,   &len);CHKERRQ(ierr);
197   ierr = PetscStrlen(ctx->bdfilename, &lenbd);CHKERRQ(ierr);
198   if (lenbd) {
199     DM bdm;
200 
201     ierr = DMPlexCreateFromFile(comm, ctx->bdfilename, PETSC_TRUE, &bdm);CHKERRQ(ierr);
202     ierr = PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_");CHKERRQ(ierr);
203     ierr = DMSetFromOptions(bdm);CHKERRQ(ierr);
204     if (ctx->scale != 0.0) {
205       Vec coordinates, coordinatesLocal;
206 
207       ierr = DMGetCoordinates(bdm, &coordinates);CHKERRQ(ierr);
208       ierr = DMGetCoordinatesLocal(bdm, &coordinatesLocal);CHKERRQ(ierr);
209       ierr = VecScale(coordinates, ctx->scale);CHKERRQ(ierr);
210       ierr = VecScale(coordinatesLocal, ctx->scale);CHKERRQ(ierr);
211     }
212     ierr = DMViewFromOptions(bdm, NULL, "-dm_view");CHKERRQ(ierr);
213     ierr = DMPlexGenerate(bdm, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);
214     ierr = DMDestroy(&bdm);CHKERRQ(ierr);
215   } else if (len) {
216     ierr = DMPlexCreateFromFile(comm, ctx->filename, PETSC_TRUE, dm);CHKERRQ(ierr);
217   } else {
218     ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);
219   }
220   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
221   ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
222   ierr = CreateBCLabel(*dm, "marker");CHKERRQ(ierr);
223   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
224   PetscFunctionReturn(0);
225 }
226 
227 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
228 {
229   PetscDS        ds;
230   const PetscInt id = 1;
231   PetscErrorCode ierr;
232 
233   PetscFunctionBeginUser;
234   ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
235   ierr = PetscDSSetJacobian(ds, 0, 0, g0_temp, NULL, NULL, g3_temp);CHKERRQ(ierr);
236   switch (ctx->solType) {
237     case SOL_QUADRATIC_LINEAR:
238       ierr = PetscDSSetResidual(ds, 0, f0_quad_lin,  f1_temp);CHKERRQ(ierr);
239       ierr = PetscDSSetExactSolution(ds, 0, mms_quad_lin, ctx);CHKERRQ(ierr);
240       ierr = PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_lin_t, ctx);CHKERRQ(ierr);
241       ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) mms_quad_lin, (void (*)(void)) mms_quad_lin_t, 1, &id, ctx);CHKERRQ(ierr);
242       break;
243     case SOL_QUADRATIC_TRIG:
244       ierr = PetscDSSetResidual(ds, 0, f0_quad_trig, f1_temp);CHKERRQ(ierr);
245       ierr = PetscDSSetExactSolution(ds, 0, mms_quad_trig, ctx);CHKERRQ(ierr);
246       ierr = PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_trig_t, ctx);CHKERRQ(ierr);
247       ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) mms_quad_trig, (void (*)(void)) mms_quad_trig_t, 1, &id, ctx);CHKERRQ(ierr);
248       break;
249     case SOL_TRIG_LINEAR:
250       ierr = PetscDSSetResidual(ds, 0, f0_trig_lin,  f1_temp);CHKERRQ(ierr);
251       ierr = PetscDSSetExactSolution(ds, 0, mms_trig_lin, ctx);CHKERRQ(ierr);
252       ierr = PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_lin_t, ctx);CHKERRQ(ierr);
253       ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) mms_trig_lin, (void (*)(void)) mms_trig_lin_t, 1, &id, ctx);CHKERRQ(ierr);
254       break;
255     default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType);
256   }
257   PetscFunctionReturn(0);
258 }
259 
260 static PetscErrorCode SetupDiscretization(DM dm, AppCtx* ctx)
261 {
262   DM             cdm = dm;
263   PetscFE        fe;
264   DMPolytopeType ct;
265   PetscBool      simplex;
266   PetscInt       dim, cStart;
267   PetscErrorCode ierr;
268 
269   PetscFunctionBeginUser;
270   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
271   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);
272   ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr);
273   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
274   /* Create finite element */
275   ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", -1, &fe);CHKERRQ(ierr);
276   ierr = PetscObjectSetName((PetscObject) fe, "temperature");CHKERRQ(ierr);
277   /* Set discretization and boundary conditions for each mesh */
278   ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr);
279   ierr = DMCreateDS(dm);CHKERRQ(ierr);
280   ierr = SetupProblem(dm, ctx);CHKERRQ(ierr);
281   while (cdm) {
282     ierr = CreateBCLabel(cdm, "marker");CHKERRQ(ierr);
283     ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr);
284     ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
285   }
286   ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
287   PetscFunctionReturn(0);
288 }
289 
290 static PetscErrorCode SetInitialConditions(TS ts, Vec u)
291 {
292   DM             dm;
293   PetscReal      t;
294   PetscErrorCode ierr;
295 
296   PetscFunctionBegin;
297   ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
298   ierr = TSGetTime(ts, &t);CHKERRQ(ierr);
299   ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr);
300   PetscFunctionReturn(0);
301 }
302 
303 int main(int argc, char **argv)
304 {
305   DM             dm;
306   TS             ts;
307   Vec            u;
308   AppCtx         ctx;
309   PetscErrorCode ierr;
310 
311   ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
312   ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr);
313   ierr = CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);CHKERRQ(ierr);
314   ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr);
315   ierr = SetupDiscretization(dm, &ctx);CHKERRQ(ierr);
316 
317   ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
318   ierr = TSSetDM(ts, dm);CHKERRQ(ierr);
319   ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);CHKERRQ(ierr);
320   ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);CHKERRQ(ierr);
321   ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);CHKERRQ(ierr);
322   ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
323   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
324   ierr = TSSetComputeInitialCondition(ts, SetInitialConditions);CHKERRQ(ierr);
325 
326   ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
327   ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr);
328   ierr = SetInitialConditions(ts, u);CHKERRQ(ierr);
329   ierr = PetscObjectSetName((PetscObject) u, "temperature");CHKERRQ(ierr);
330   ierr = TSSolve(ts, u);CHKERRQ(ierr);
331   ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr);
332 
333   ierr = VecDestroy(&u);CHKERRQ(ierr);
334   ierr = TSDestroy(&ts);CHKERRQ(ierr);
335   ierr = DMDestroy(&dm);CHKERRQ(ierr);
336   ierr = PetscFinalize();
337   return ierr;
338 }
339 
340 /*TEST
341 
342   test:
343     suffix: 2d_p1
344     requires: triangle
345     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
346           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
347   test:
348     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
349     suffix: 2d_p1_sconv
350     requires: triangle
351     args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
352           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
353   test:
354     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
355     suffix: 2d_p1_tconv
356     requires: triangle
357     args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
358           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
359   test:
360     suffix: 2d_p2
361     requires: triangle
362     args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -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     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
366     suffix: 2d_p2_sconv
367     requires: triangle
368     args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
369           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
370   test:
371     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
372     suffix: 2d_p2_tconv
373     requires: triangle
374     args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
375           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
376   test:
377     suffix: 2d_q1
378     args: -sol_type quadratic_linear -dm_plex_box_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
379           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
380   test:
381     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
382     suffix: 2d_q1_sconv
383     args: -sol_type quadratic_linear -dm_plex_box_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
384           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
385   test:
386     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
387     suffix: 2d_q1_tconv
388     args: -sol_type quadratic_trig -dm_plex_box_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
389           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
390   test:
391     suffix: 2d_q2
392     args: -sol_type quadratic_linear -dm_plex_box_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -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     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
396     suffix: 2d_q2_sconv
397     args: -sol_type trig_linear -dm_plex_box_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
398           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
399   test:
400     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
401     suffix: 2d_q2_tconv
402     args: -sol_type quadratic_trig -dm_plex_box_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
403           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
404 
405   test:
406     suffix: 3d_p1
407     requires: ctetgen
408     args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -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: [1.9]
412     suffix: 3d_p1_sconv
413     requires: ctetgen
414     args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
415           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
416   test:
417     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
418     suffix: 3d_p1_tconv
419     requires: ctetgen
420     args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
421           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
422   test:
423     suffix: 3d_p2
424     requires: ctetgen
425     args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
426           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
427   test:
428     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
429     suffix: 3d_p2_sconv
430     requires: ctetgen
431     args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
432           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
433   test:
434     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
435     suffix: 3d_p2_tconv
436     requires: ctetgen
437     args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
438           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
439   test:
440     suffix: 3d_q1
441     args: -sol_type quadratic_linear -dm_plex_box_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \
442           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
443   test:
444     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9]
445     suffix: 3d_q1_sconv
446     args: -sol_type quadratic_linear -dm_plex_box_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
447           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu
448   test:
449     # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2]
450     suffix: 3d_q1_tconv
451     args: -sol_type quadratic_trig -dm_plex_box_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \
452           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
453   test:
454     suffix: 3d_q2
455     args: -sol_type quadratic_linear -dm_plex_box_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \
456           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
457   test:
458     # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9]
459     suffix: 3d_q2_sconv
460     args: -sol_type trig_linear -dm_plex_box_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
461           -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu
462   test:
463     # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0]
464     suffix: 3d_q2_tconv
465     args: -sol_type quadratic_trig -dm_plex_box_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \
466           -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
467 
468   test:
469     # 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
470     suffix: egads_sphere
471     requires: egads ctetgen
472     args: -sol_type quadratic_linear -bd_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -scale 40 \
473           -temp_petscspace_degree 2 -dmts_check .0001 \
474           -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu
475 
476 TEST*/
477