xref: /petsc/src/dm/impls/plex/tests/ex9.c (revision d5b43468fb8780a8feea140ccd6fa3e6a50411cc)
1 static char help[] = "Performance tests for DMPlex query operations\n\n";
2 
3 #include <petscdmplex.h>
4 
5 typedef struct {
6   PetscInt  dim;             /* The topological mesh dimension */
7   PetscBool cellSimplex;     /* Flag for simplices */
8   PetscBool spectral;        /* Flag for spectral element layout */
9   PetscBool interpolate;     /* Flag for mesh interpolation */
10   PetscReal refinementLimit; /* Maximum volume of a refined cell */
11   PetscInt  numFields;       /* The number of section fields */
12   PetscInt *numComponents;   /* The number of field components */
13   PetscInt *numDof;          /* The dof signature for the section */
14   PetscBool reuseArray;      /* Pass in user allocated array to VecGetClosure() */
15   /* Test data */
16   PetscBool errors;            /* Treat failures as errors */
17   PetscInt  iterations;        /* The number of iterations for a query */
18   PetscReal maxConeTime;       /* Max time per run for DMPlexGetCone() */
19   PetscReal maxClosureTime;    /* Max time per run for DMPlexGetTransitiveClosure() */
20   PetscReal maxVecClosureTime; /* Max time per run for DMPlexVecGetClosure() */
21   PetscBool printTimes;        /* Print total times, do not check limits */
22 } AppCtx;
23 
24 static PetscErrorCode ProcessOptions(AppCtx *options)
25 {
26   PetscInt  len;
27   PetscBool flg;
28 
29   PetscFunctionBegin;
30   options->dim               = 2;
31   options->cellSimplex       = PETSC_TRUE;
32   options->spectral          = PETSC_FALSE;
33   options->interpolate       = PETSC_FALSE;
34   options->refinementLimit   = 0.0;
35   options->numFields         = 0;
36   options->numComponents     = NULL;
37   options->numDof            = NULL;
38   options->reuseArray        = PETSC_FALSE;
39   options->errors            = PETSC_FALSE;
40   options->iterations        = 1;
41   options->maxConeTime       = 0.0;
42   options->maxClosureTime    = 0.0;
43   options->maxVecClosureTime = 0.0;
44   options->printTimes        = PETSC_FALSE;
45 
46   PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
47   PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex9.c", options->dim, &options->dim, NULL, 1, 3));
48   PetscCall(PetscOptionsBool("-cellSimplex", "Flag for simplices", "ex9.c", options->cellSimplex, &options->cellSimplex, NULL));
49   PetscCall(PetscOptionsBool("-spectral", "Flag for spectral element layout", "ex9.c", options->spectral, &options->spectral, NULL));
50   PetscCall(PetscOptionsBool("-interpolate", "Flag for mesh interpolation", "ex9.c", options->interpolate, &options->interpolate, NULL));
51   PetscCall(PetscOptionsReal("-refinement_limit", "The maximum volume of a refined cell", "ex9.c", options->refinementLimit, &options->refinementLimit, NULL));
52   PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of section fields", "ex9.c", options->numFields, &options->numFields, NULL, 0));
53   if (options->numFields) {
54     len = options->numFields;
55     PetscCall(PetscMalloc1(len, &options->numComponents));
56     PetscCall(PetscOptionsIntArray("-num_components", "The number of components per field", "ex9.c", options->numComponents, &len, &flg));
57     PetscCheck(!flg || !(len != options->numFields), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of components array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, options->numFields);
58   }
59   len = (options->dim + 1) * PetscMax(1, options->numFields);
60   PetscCall(PetscMalloc1(len, &options->numDof));
61   PetscCall(PetscOptionsIntArray("-num_dof", "The dof signature for the section", "ex9.c", options->numDof, &len, &flg));
62   PetscCheck(!flg || len == (options->dim + 1) * PetscMax(1, options->numFields), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of dof array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, (options->dim + 1) * PetscMax(1, options->numFields));
63 
64   /* We are specifying the scalar dof, so augment it for multiple components */
65   {
66     PetscInt f, d;
67 
68     for (f = 0; f < options->numFields; ++f) {
69       for (d = 0; d <= options->dim; ++d) options->numDof[f * (options->dim + 1) + d] *= options->numComponents[f];
70     }
71   }
72 
73   PetscCall(PetscOptionsBool("-reuse_array", "Pass in user allocated array to VecGetClosure()", "ex9.c", options->reuseArray, &options->reuseArray, NULL));
74   PetscCall(PetscOptionsBool("-errors", "Treat failures as errors", "ex9.c", options->errors, &options->errors, NULL));
75   PetscCall(PetscOptionsBoundedInt("-iterations", "The number of iterations for a query", "ex9.c", options->iterations, &options->iterations, NULL, 0));
76   PetscCall(PetscOptionsReal("-max_cone_time", "The maximum time per run for DMPlexGetCone()", "ex9.c", options->maxConeTime, &options->maxConeTime, NULL));
77   PetscCall(PetscOptionsReal("-max_closure_time", "The maximum time per run for DMPlexGetTransitiveClosure()", "ex9.c", options->maxClosureTime, &options->maxClosureTime, NULL));
78   PetscCall(PetscOptionsReal("-max_vec_closure_time", "The maximum time per run for DMPlexVecGetClosure()", "ex9.c", options->maxVecClosureTime, &options->maxVecClosureTime, NULL));
79   PetscCall(PetscOptionsBool("-print_times", "Print total times, do not check limits", "ex9.c", options->printTimes, &options->printTimes, NULL));
80   PetscOptionsEnd();
81   PetscFunctionReturn(0);
82 }
83 
84 static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, DM *newdm)
85 {
86   DM          dm;
87   PetscInt    numPoints[2]        = {4, 2};
88   PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
89   PetscInt    cones[6]            = {2, 3, 4, 5, 4, 3};
90   PetscInt    coneOrientations[6] = {0, 0, 0, 0, 0, 0};
91   PetscScalar vertexCoords[8]     = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
92   PetscInt    markerPoints[8]     = {2, 1, 3, 1, 4, 1, 5, 1};
93   PetscInt    dim = 2, depth = 1, p;
94 
95   PetscFunctionBegin;
96   PetscCall(DMCreate(comm, &dm));
97   PetscCall(PetscObjectSetName((PetscObject)dm, "triangular"));
98   PetscCall(DMSetType(dm, DMPLEX));
99   PetscCall(DMSetDimension(dm, dim));
100   PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
101   for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
102   *newdm = dm;
103   PetscFunctionReturn(0);
104 }
105 
106 static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, DM *newdm)
107 {
108   DM          dm;
109   PetscInt    numPoints[2]        = {5, 2};
110   PetscInt    coneSize[23]        = {4, 4, 0, 0, 0, 0, 0};
111   PetscInt    cones[8]            = {2, 4, 3, 5, 3, 4, 6, 5};
112   PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
113   PetscScalar vertexCoords[15]    = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5};
114   PetscInt    markerPoints[10]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
115   PetscInt    dim = 3, depth = 1, p;
116 
117   PetscFunctionBegin;
118   PetscCall(DMCreate(comm, &dm));
119   PetscCall(PetscObjectSetName((PetscObject)dm, "tetrahedral"));
120   PetscCall(DMSetType(dm, DMPLEX));
121   PetscCall(DMSetDimension(dm, dim));
122   PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
123   for (p = 0; p < 5; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
124   *newdm = dm;
125   PetscFunctionReturn(0);
126 }
127 
128 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, DM *newdm)
129 {
130   DM          dm;
131   PetscInt    numPoints[2]        = {6, 2};
132   PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
133   PetscInt    cones[8]            = {2, 3, 4, 5, 3, 6, 7, 4};
134   PetscInt    coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
135   PetscScalar vertexCoords[12]    = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0, 0.5, 1.0};
136   PetscInt    markerPoints[12]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
137   PetscInt    dim = 2, depth = 1, p;
138 
139   PetscFunctionBegin;
140   PetscCall(DMCreate(comm, &dm));
141   PetscCall(PetscObjectSetName((PetscObject)dm, "quadrilateral"));
142   PetscCall(DMSetType(dm, DMPLEX));
143   PetscCall(DMSetDimension(dm, dim));
144   PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
145   for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
146   *newdm = dm;
147   PetscFunctionReturn(0);
148 }
149 
150 static PetscErrorCode CreateHex_3D(MPI_Comm comm, DM *newdm)
151 {
152   DM          dm;
153   PetscInt    numPoints[2]         = {12, 2};
154   PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
155   PetscInt    cones[16]            = {2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8};
156   PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
157   PetscScalar vertexCoords[36]     = {-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -0.5, 1.0, 0.0, -0.5, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, -0.5, 1.0, 1.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0};
158   PetscInt    markerPoints[24]     = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1, 8, 1, 9, 1, 10, 1, 11, 1, 12, 1, 13, 1};
159   PetscInt    dim = 3, depth = 1, p;
160 
161   PetscFunctionBegin;
162   PetscCall(DMCreate(comm, &dm));
163   PetscCall(PetscObjectSetName((PetscObject)dm, "hexahedral"));
164   PetscCall(DMSetType(dm, DMPLEX));
165   PetscCall(DMSetDimension(dm, dim));
166   PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords));
167   for (p = 0; p < 12; ++p) PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p * 2], markerPoints[p * 2 + 1]));
168   *newdm = dm;
169   PetscFunctionReturn(0);
170 }
171 
172 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *newdm)
173 {
174   PetscInt  dim         = user->dim;
175   PetscBool cellSimplex = user->cellSimplex;
176 
177   PetscFunctionBegin;
178   switch (dim) {
179   case 2:
180     if (cellSimplex) {
181       PetscCall(CreateSimplex_2D(comm, newdm));
182     } else {
183       PetscCall(CreateQuad_2D(comm, newdm));
184     }
185     break;
186   case 3:
187     if (cellSimplex) {
188       PetscCall(CreateSimplex_3D(comm, newdm));
189     } else {
190       PetscCall(CreateHex_3D(comm, newdm));
191     }
192     break;
193   default:
194     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %" PetscInt_FMT, dim);
195   }
196   if (user->refinementLimit > 0.0) {
197     DM          rdm;
198     const char *name;
199 
200     PetscCall(DMPlexSetRefinementUniform(*newdm, PETSC_FALSE));
201     PetscCall(DMPlexSetRefinementLimit(*newdm, user->refinementLimit));
202     PetscCall(DMRefine(*newdm, PETSC_COMM_SELF, &rdm));
203     PetscCall(PetscObjectGetName((PetscObject)*newdm, &name));
204     PetscCall(PetscObjectSetName((PetscObject)rdm, name));
205     PetscCall(DMDestroy(newdm));
206     *newdm = rdm;
207   }
208   if (user->interpolate) {
209     DM idm;
210 
211     PetscCall(DMPlexInterpolate(*newdm, &idm));
212     PetscCall(DMDestroy(newdm));
213     *newdm = idm;
214   }
215   PetscCall(DMSetFromOptions(*newdm));
216   PetscFunctionReturn(0);
217 }
218 
219 static PetscErrorCode TestCone(DM dm, AppCtx *user)
220 {
221   PetscInt           numRuns, cStart, cEnd, c, i;
222   PetscReal          maxTimePerRun = user->maxConeTime;
223   PetscLogStage      stage;
224   PetscLogEvent      event;
225   PetscEventPerfInfo eventInfo;
226   MPI_Comm           comm;
227   PetscMPIInt        rank;
228 
229   PetscFunctionBegin;
230   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
231   PetscCallMPI(MPI_Comm_rank(comm, &rank));
232   PetscCall(PetscLogStageRegister("DMPlex Cone Test", &stage));
233   PetscCall(PetscLogEventRegister("Cone", PETSC_OBJECT_CLASSID, &event));
234   PetscCall(PetscLogStagePush(stage));
235   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
236   PetscCall(PetscLogEventBegin(event, 0, 0, 0, 0));
237   for (i = 0; i < user->iterations; ++i) {
238     for (c = cStart; c < cEnd; ++c) {
239       const PetscInt *cone;
240 
241       PetscCall(DMPlexGetCone(dm, c, &cone));
242     }
243   }
244   PetscCall(PetscLogEventEnd(event, 0, 0, 0, 0));
245   PetscCall(PetscLogStagePop());
246 
247   PetscCall(PetscLogEventGetPerfInfo(stage, event, &eventInfo));
248   numRuns = (cEnd - cStart) * user->iterations;
249   PetscCheck(eventInfo.count == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be 1", eventInfo.count);
250   PetscCheck((PetscInt)eventInfo.flops == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %" PetscInt_FMT " should be 0", (PetscInt)eventInfo.flops);
251   if (user->printTimes) {
252     PetscCall(PetscSynchronizedPrintf(comm, "[%d] Cones: %" PetscInt_FMT " Total time: %.3es Average time per cone: %.3es\n", rank, numRuns, eventInfo.time, eventInfo.time / numRuns));
253     PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
254   } else if (eventInfo.time > maxTimePerRun * numRuns) {
255     PetscCall(PetscSynchronizedPrintf(comm, "[%d] Cones: %" PetscInt_FMT " Average time per cone: %gs standard: %gs\n", rank, numRuns, (double)(eventInfo.time / numRuns), (double)maxTimePerRun));
256     PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
257     PetscCheck(!user->errors, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for cone %g > standard %g", (double)(eventInfo.time / numRuns), (double)maxTimePerRun);
258   }
259   PetscFunctionReturn(0);
260 }
261 
262 static PetscErrorCode TestTransitiveClosure(DM dm, AppCtx *user)
263 {
264   PetscInt           numRuns, cStart, cEnd, c, i;
265   PetscReal          maxTimePerRun = user->maxClosureTime;
266   PetscLogStage      stage;
267   PetscLogEvent      event;
268   PetscEventPerfInfo eventInfo;
269   MPI_Comm           comm;
270   PetscMPIInt        rank;
271 
272   PetscFunctionBegin;
273   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
274   PetscCallMPI(MPI_Comm_rank(comm, &rank));
275   PetscCall(PetscLogStageRegister("DMPlex Transitive Closure Test", &stage));
276   PetscCall(PetscLogEventRegister("TransitiveClosure", PETSC_OBJECT_CLASSID, &event));
277   PetscCall(PetscLogStagePush(stage));
278   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
279   PetscCall(PetscLogEventBegin(event, 0, 0, 0, 0));
280   for (i = 0; i < user->iterations; ++i) {
281     for (c = cStart; c < cEnd; ++c) {
282       PetscInt *closure = NULL;
283       PetscInt  closureSize;
284 
285       PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
286       PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
287     }
288   }
289   PetscCall(PetscLogEventEnd(event, 0, 0, 0, 0));
290   PetscCall(PetscLogStagePop());
291 
292   PetscCall(PetscLogEventGetPerfInfo(stage, event, &eventInfo));
293   numRuns = (cEnd - cStart) * user->iterations;
294   PetscCheck(eventInfo.count == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be 1", eventInfo.count);
295   PetscCheck((PetscInt)eventInfo.flops == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %" PetscInt_FMT " should be 0", (PetscInt)eventInfo.flops);
296   if (user->printTimes) {
297     PetscCall(PetscSynchronizedPrintf(comm, "[%d] Closures: %" PetscInt_FMT " Total time: %.3es Average time per cone: %.3es\n", rank, numRuns, eventInfo.time, eventInfo.time / numRuns));
298     PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
299   } else if (eventInfo.time > maxTimePerRun * numRuns) {
300     PetscCall(PetscSynchronizedPrintf(comm, "[%d] Closures: %" PetscInt_FMT " Average time per cone: %gs standard: %gs\n", rank, numRuns, (double)(eventInfo.time / numRuns), (double)maxTimePerRun));
301     PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
302     PetscCheck(!user->errors, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for closure %g > standard %g", (double)(eventInfo.time / numRuns), (double)maxTimePerRun);
303   }
304   PetscFunctionReturn(0);
305 }
306 
307 static PetscErrorCode TestVecClosure(DM dm, PetscBool useIndex, PetscBool useSpectral, AppCtx *user)
308 {
309   PetscSection       s;
310   Vec                v;
311   PetscInt           numRuns, cStart, cEnd, c, i;
312   PetscScalar        tmpArray[64];
313   PetscScalar       *userArray     = user->reuseArray ? tmpArray : NULL;
314   PetscReal          maxTimePerRun = user->maxVecClosureTime;
315   PetscLogStage      stage;
316   PetscLogEvent      event;
317   PetscEventPerfInfo eventInfo;
318   MPI_Comm           comm;
319   PetscMPIInt        rank;
320 
321   PetscFunctionBegin;
322   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
323   PetscCallMPI(MPI_Comm_rank(comm, &rank));
324   if (useIndex) {
325     if (useSpectral) {
326       PetscCall(PetscLogStageRegister("DMPlex Vector Closure with Index Test", &stage));
327       PetscCall(PetscLogEventRegister("VecClosureInd", PETSC_OBJECT_CLASSID, &event));
328     } else {
329       PetscCall(PetscLogStageRegister("DMPlex Vector Spectral Closure with Index Test", &stage));
330       PetscCall(PetscLogEventRegister("VecClosureSpecInd", PETSC_OBJECT_CLASSID, &event));
331     }
332   } else {
333     if (useSpectral) {
334       PetscCall(PetscLogStageRegister("DMPlex Vector Spectral Closure Test", &stage));
335       PetscCall(PetscLogEventRegister("VecClosureSpec", PETSC_OBJECT_CLASSID, &event));
336     } else {
337       PetscCall(PetscLogStageRegister("DMPlex Vector Closure Test", &stage));
338       PetscCall(PetscLogEventRegister("VecClosure", PETSC_OBJECT_CLASSID, &event));
339     }
340   }
341   PetscCall(PetscLogStagePush(stage));
342   PetscCall(DMSetNumFields(dm, user->numFields));
343   PetscCall(DMPlexCreateSection(dm, NULL, user->numComponents, user->numDof, 0, NULL, NULL, NULL, NULL, &s));
344   PetscCall(DMSetLocalSection(dm, s));
345   if (useIndex) PetscCall(DMPlexCreateClosureIndex(dm, s));
346   if (useSpectral) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, s));
347   PetscCall(PetscSectionDestroy(&s));
348   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
349   PetscCall(DMGetLocalVector(dm, &v));
350   PetscCall(PetscLogEventBegin(event, 0, 0, 0, 0));
351   for (i = 0; i < user->iterations; ++i) {
352     for (c = cStart; c < cEnd; ++c) {
353       PetscScalar *closure     = userArray;
354       PetscInt     closureSize = 64;
355 
356       PetscCall(DMPlexVecGetClosure(dm, s, v, c, &closureSize, &closure));
357       if (!user->reuseArray) PetscCall(DMPlexVecRestoreClosure(dm, s, v, c, &closureSize, &closure));
358     }
359   }
360   PetscCall(PetscLogEventEnd(event, 0, 0, 0, 0));
361   PetscCall(DMRestoreLocalVector(dm, &v));
362   PetscCall(PetscLogStagePop());
363 
364   PetscCall(PetscLogEventGetPerfInfo(stage, event, &eventInfo));
365   numRuns = (cEnd - cStart) * user->iterations;
366   PetscCheck(eventInfo.count == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be 1", eventInfo.count);
367   PetscCheck((PetscInt)eventInfo.flops == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %" PetscInt_FMT " should be 0", (PetscInt)eventInfo.flops);
368   if (user->printTimes || eventInfo.time > maxTimePerRun * numRuns) {
369     const char *title          = "VecClosures";
370     const char *titleIndex     = "VecClosures with Index";
371     const char *titleSpec      = "VecClosures Spectral";
372     const char *titleSpecIndex = "VecClosures Spectral with Index";
373 
374     if (user->printTimes) {
375       PetscCall(PetscSynchronizedPrintf(comm, "[%d] %s: %" PetscInt_FMT " Total time: %.3es Average time per vector closure: %.3es\n", rank, useIndex ? (useSpectral ? titleSpecIndex : titleIndex) : (useSpectral ? titleSpec : title), numRuns,
376                                         eventInfo.time, eventInfo.time / numRuns));
377       PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
378     } else {
379       PetscCall(
380         PetscSynchronizedPrintf(comm, "[%d] %s: %" PetscInt_FMT " Average time per vector closure: %gs standard: %gs\n", rank, useIndex ? (useSpectral ? titleSpecIndex : titleIndex) : (useSpectral ? titleSpec : title), numRuns, (double)(eventInfo.time / numRuns), (double)maxTimePerRun));
381       PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
382       PetscCheck(!user->errors, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for vector closure %g > standard %g", (double)(eventInfo.time / numRuns), (double)maxTimePerRun);
383     }
384   }
385   PetscFunctionReturn(0);
386 }
387 
388 static PetscErrorCode CleanupContext(AppCtx *user)
389 {
390   PetscFunctionBegin;
391   PetscCall(PetscFree(user->numComponents));
392   PetscCall(PetscFree(user->numDof));
393   PetscFunctionReturn(0);
394 }
395 
396 int main(int argc, char **argv)
397 {
398   DM     dm;
399   AppCtx user;
400 
401   PetscFunctionBeginUser;
402   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
403   PetscCall(ProcessOptions(&user));
404   PetscCall(PetscLogDefaultBegin());
405   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
406   PetscCall(TestCone(dm, &user));
407   PetscCall(TestTransitiveClosure(dm, &user));
408   PetscCall(TestVecClosure(dm, PETSC_FALSE, PETSC_FALSE, &user));
409   PetscCall(TestVecClosure(dm, PETSC_TRUE, PETSC_FALSE, &user));
410   if (!user.cellSimplex && user.spectral) {
411     PetscCall(TestVecClosure(dm, PETSC_FALSE, PETSC_TRUE, &user));
412     PetscCall(TestVecClosure(dm, PETSC_TRUE, PETSC_TRUE, &user));
413   }
414   PetscCall(DMDestroy(&dm));
415   PetscCall(CleanupContext(&user));
416   PetscCall(PetscFinalize());
417   return 0;
418 }
419 
420 /*TEST
421 
422   build:
423     requires: defined(PETSC_USE_LOG)
424 
425   # 2D Simplex P_1 scalar tests
426   testset:
427     args: -num_dof 1,0,0 -iterations 2 -print_times
428     test:
429       suffix: correctness_0
430     test:
431       suffix: correctness_1
432       args: -interpolate -dm_refine 2
433     test:
434       suffix: correctness_2
435       requires: triangle
436       args: -interpolate -refinement_limit 1.0e-5
437   test:
438     suffix: 0
439     TODO: Only for performance testing
440     args: -num_dof 1,0,0 -iterations 10000 -max_cone_time 1.1e-8 -max_closure_time 1.3e-7 -max_vec_closure_time 3.6e-7
441   test:
442     suffix: 1
443     requires: triangle
444     TODO: Only for performance testing
445     args: -refinement_limit 1.0e-5 -num_dof 1,0,0 -iterations 2 -max_cone_time 2.1e-8 -max_closure_time 1.5e-7 -max_vec_closure_time 3.6e-7
446   test:
447     suffix: 2
448     TODO: Only for performance testing
449     args: -num_fields 1 -num_components 1 -num_dof 1,0,0 -iterations 10000 -max_cone_time 1.1e-8 -max_closure_time 1.3e-7 -max_vec_closure_time 4.5e-7
450   test:
451     suffix: 3
452     requires: triangle
453     TODO: Only for performance testing
454     args: -refinement_limit 1.0e-5 -num_fields 1 -num_components 1 -num_dof 1,0,0 -iterations 2 -max_cone_time 2.1e-8 -max_closure_time 1.5e-7 -max_vec_closure_time 4.7e-7
455   test:
456     suffix: 4
457     TODO: Only for performance testing
458     args: -interpolate -num_dof 1,0,0 -iterations 10000 -max_cone_time 1.1e-8 -max_closure_time 6.5e-7 -max_vec_closure_time 1.0e-6
459   test:
460     suffix: 5
461     requires: triangle
462     TODO: Only for performance testing
463     args: -interpolate -refinement_limit 1.0e-4 -num_dof 1,0,0 -iterations 2 -max_cone_time 2.1e-8 -max_closure_time 6.5e-7 -max_vec_closure_time 1.0e-6
464   test:
465     suffix: 6
466     TODO: Only for performance testing
467     args: -interpolate -num_fields 1 -num_components 1 -num_dof 1,0,0 -iterations 10000 -max_cone_time 1.1e-8 -max_closure_time 6.5e-7 -max_vec_closure_time 1.1e-6
468   test:
469     suffix: 7
470     requires: triangle
471     TODO: Only for performance testing
472     args: -interpolate -refinement_limit 1.0e-4 -num_fields 1 -num_components 1 -num_dof 1,0,0 -iterations 2 -max_cone_time 2.1e-8 -max_closure_time 6.5e-7 -max_vec_closure_time 1.2e-6
473 
474   # 2D Simplex P_1 vector tests
475   # 2D Simplex P_2 scalar tests
476   # 2D Simplex P_2 vector tests
477   # 2D Simplex P_2/P_1 vector/scalar tests
478   # 2D Quad P_1 scalar tests
479   # 2D Quad P_1 vector tests
480   # 2D Quad P_2 scalar tests
481   # 2D Quad P_2 vector tests
482   # 3D Simplex P_1 scalar tests
483   # 3D Simplex P_1 vector tests
484   # 3D Simplex P_2 scalar tests
485   # 3D Simplex P_2 vector tests
486   # 3D Hex P_1 scalar tests
487   # 3D Hex P_1 vector tests
488   # 3D Hex P_2 scalar tests
489   # 3D Hex P_2 vector tests
490 
491 TEST*/
492