Actual source code: ex26.c
1: static char help[] = "Test FEM layout with DM and ExodusII storage\n\n";
3: /*
4: In order to see the vectors which are being tested, use
6: -ua_vec_view -s_vec_view
7: */
9: #include <petsc.h>
10: #include <exodusII.h>
12: #include <petsc/private/dmpleximpl.h>
14: int main(int argc, char **argv)
15: {
16: DM dm, pdm, dmU, dmA, dmS, dmUA, dmUA2, *dmList;
17: Vec X, U, A, S, UA, UA2;
18: IS isU, isA, isS, isUA;
19: PetscSection section;
20: const PetscInt fieldU = 0;
21: const PetscInt fieldA = 2;
22: const PetscInt fieldS = 1;
23: const PetscInt fieldUA[2] = {0, 2};
24: char ifilename[PETSC_MAX_PATH_LEN], ofilename[PETSC_MAX_PATH_LEN];
25: int exoid = -1;
26: IS csIS;
27: const PetscInt *csID;
28: PetscInt *pStartDepth, *pEndDepth;
29: PetscInt order = 1;
30: PetscInt sdim, d, pStart, pEnd, p, numCS, set;
31: PetscMPIInt rank, size;
32: PetscViewer viewer;
34: PetscFunctionBeginUser;
35: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
36: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
37: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
38: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FEM Layout Options", "ex26");
39: PetscCall(PetscOptionsString("-i", "Filename to read", "ex26", ifilename, ifilename, sizeof(ifilename), NULL));
40: PetscCall(PetscOptionsString("-o", "Filename to write", "ex26", ofilename, ofilename, sizeof(ofilename), NULL));
41: PetscCall(PetscOptionsBoundedInt("-order", "FEM polynomial order", "ex26", order, &order, NULL, 1));
42: PetscOptionsEnd();
43: PetscCheck((order >= 1) && (order <= 2), PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported polynomial order %" PetscInt_FMT " not in [1, 2]", order);
45: /* Read the mesh from a file in any supported format */
46: PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, NULL, PETSC_TRUE, &dm));
47: PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE));
48: PetscCall(DMSetFromOptions(dm));
49: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
50: PetscCall(DMGetDimension(dm, &sdim));
52: /* Create the exodus result file */
53: {
54: PetscInt numstep = 3, step;
55: char *nodalVarName[4];
56: char *zonalVarName[6];
57: int *truthtable;
58: PetscInt numNodalVar, numZonalVar, i;
60: /* enable exodus debugging information */
61: ex_opts(EX_VERBOSE | EX_DEBUG);
62: /* Create the exodus file */
63: PetscCall(PetscViewerExodusIIOpen(PETSC_COMM_WORLD, ofilename, FILE_MODE_WRITE, &viewer));
64: /* The long way would be */
65: /*
66: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD,&viewer));
67: PetscCall(PetscViewerSetType(viewer,PETSCVIEWEREXODUSII));
68: PetscCall(PetscViewerFileSetMode(viewer,FILE_MODE_APPEND));
69: PetscCall(PetscViewerFileSetName(viewer,ofilename));
70: */
71: /* set the mesh order */
72: PetscCall(PetscViewerExodusIISetOrder(viewer, order));
73: PetscCall(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD));
74: /*
75: Notice how the exodus file is actually NOT open at this point (exoid is -1)
76: Since we are overwriting the file (mode is FILE_MODE_WRITE), we are going to have to
77: write the geometry (the DM), which can only be done on a brand new file.
78: */
80: /* Save the geometry to the file, erasing all previous content */
81: PetscCall(DMView(dm, viewer));
82: PetscCall(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD));
83: /*
84: Note how the exodus file is now open
85: */
87: /* "Format" the exodus result file, i.e. allocate space for nodal and zonal variables */
88: switch (sdim) {
89: case 2:
90: numNodalVar = 3;
91: nodalVarName[0] = (char *)"U_x";
92: nodalVarName[1] = (char *)"U_y";
93: nodalVarName[2] = (char *)"Alpha";
94: numZonalVar = 3;
95: zonalVarName[0] = (char *)"Sigma_11";
96: zonalVarName[1] = (char *)"Sigma_22";
97: zonalVarName[2] = (char *)"Sigma_12";
98: break;
99: case 3:
100: numNodalVar = 4;
101: nodalVarName[0] = (char *)"U_x";
102: nodalVarName[1] = (char *)"U_y";
103: nodalVarName[2] = (char *)"U_z";
104: nodalVarName[3] = (char *)"Alpha";
105: numZonalVar = 6;
106: zonalVarName[0] = (char *)"Sigma_11";
107: zonalVarName[1] = (char *)"Sigma_22";
108: zonalVarName[2] = (char *)"Sigma_33";
109: zonalVarName[3] = (char *)"Sigma_23";
110: zonalVarName[4] = (char *)"Sigma_13";
111: zonalVarName[5] = (char *)"Sigma_12";
112: break;
113: default:
114: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %" PetscInt_FMT, sdim);
115: }
116: PetscCall(PetscViewerExodusIIGetId(viewer, &exoid));
117: PetscCallExternal(ex_put_variable_param, exoid, EX_ELEM_BLOCK, numZonalVar);
118: PetscCallExternal(ex_put_variable_names, exoid, EX_ELEM_BLOCK, numZonalVar, zonalVarName);
119: PetscCallExternal(ex_put_variable_param, exoid, EX_NODAL, numNodalVar);
120: PetscCallExternal(ex_put_variable_names, exoid, EX_NODAL, numNodalVar, nodalVarName);
121: numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
123: /*
124: An exodusII truth table specifies which fields are saved at which time step
125: It speeds up I/O but reserving space for fieldsin the file ahead of time.
126: */
127: PetscCall(PetscMalloc1(numZonalVar * numCS, &truthtable));
128: for (i = 0; i < numZonalVar * numCS; ++i) truthtable[i] = 1;
129: PetscCallExternal(ex_put_truth_table, exoid, EX_ELEM_BLOCK, numCS, numZonalVar, truthtable);
130: PetscCall(PetscFree(truthtable));
132: /* Writing time step information in the file. Note that this is currently broken in the exodus library for netcdf4 (HDF5-based) files */
133: for (step = 0; step < numstep; ++step) {
134: PetscReal time = step;
135: PetscCallExternal(ex_put_time, exoid, step + 1, &time);
136: }
137: }
139: /* Create the main section containing all fields */
140: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
141: PetscCall(PetscSectionSetNumFields(section, 3));
142: PetscCall(PetscSectionSetFieldName(section, fieldU, "U"));
143: PetscCall(PetscSectionSetFieldName(section, fieldA, "Alpha"));
144: PetscCall(PetscSectionSetFieldName(section, fieldS, "Sigma"));
145: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
146: PetscCall(PetscSectionSetChart(section, pStart, pEnd));
147: PetscCall(PetscMalloc2(sdim + 1, &pStartDepth, sdim + 1, &pEndDepth));
148: for (d = 0; d <= sdim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStartDepth[d], &pEndDepth[d]));
149: /* Vector field U, Scalar field Alpha, Tensor field Sigma */
150: PetscCall(PetscSectionSetFieldComponents(section, fieldU, sdim));
151: PetscCall(PetscSectionSetFieldComponents(section, fieldA, 1));
152: PetscCall(PetscSectionSetFieldComponents(section, fieldS, sdim * (sdim + 1) / 2));
154: /* Going through cell sets then cells, and setting up storage for the sections */
155: PetscCall(DMGetLabelSize(dm, "Cell Sets", &numCS));
156: PetscCall(DMGetLabelIdIS(dm, "Cell Sets", &csIS));
157: if (csIS) PetscCall(ISGetIndices(csIS, &csID));
158: for (set = 0; set < numCS; set++) {
159: IS cellIS;
160: const PetscInt *cellID;
161: PetscInt numCells, cell, closureSize, *closureA = NULL;
163: PetscCall(DMGetStratumSize(dm, "Cell Sets", csID[set], &numCells));
164: PetscCall(DMGetStratumIS(dm, "Cell Sets", csID[set], &cellIS));
165: if (numCells > 0) {
166: /* dof layout ordered by increasing height in the DAG: cell, face, edge, vertex */
167: PetscInt dofUP1Tri[] = {2, 0, 0};
168: PetscInt dofAP1Tri[] = {1, 0, 0};
169: PetscInt dofUP2Tri[] = {2, 2, 0};
170: PetscInt dofAP2Tri[] = {1, 1, 0};
171: PetscInt dofUP1Quad[] = {2, 0, 0};
172: PetscInt dofAP1Quad[] = {1, 0, 0};
173: PetscInt dofUP2Quad[] = {2, 2, 2};
174: PetscInt dofAP2Quad[] = {1, 1, 1};
175: PetscInt dofS2D[] = {0, 0, 3};
176: PetscInt dofUP1Tet[] = {3, 0, 0, 0};
177: PetscInt dofAP1Tet[] = {1, 0, 0, 0};
178: PetscInt dofUP2Tet[] = {3, 3, 0, 0};
179: PetscInt dofAP2Tet[] = {1, 1, 0, 0};
180: PetscInt dofUP1Hex[] = {3, 0, 0, 0};
181: PetscInt dofAP1Hex[] = {1, 0, 0, 0};
182: PetscInt dofUP2Hex[] = {3, 3, 3, 3};
183: PetscInt dofAP2Hex[] = {1, 1, 1, 1};
184: PetscInt dofS3D[] = {0, 0, 0, 6};
185: PetscInt *dofU, *dofA, *dofS;
187: switch (sdim) {
188: case 2:
189: dofS = dofS2D;
190: break;
191: case 3:
192: dofS = dofS3D;
193: break;
194: default:
195: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %" PetscInt_FMT, sdim);
196: }
198: /* Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes
199: It will not be enough to identify more exotic elements like pyramid or prisms... */
200: PetscCall(ISGetIndices(cellIS, &cellID));
201: PetscCall(DMPlexGetTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA));
202: switch (closureSize) {
203: case 7: /* Tri */
204: if (order == 1) {
205: dofU = dofUP1Tri;
206: dofA = dofAP1Tri;
207: } else {
208: dofU = dofUP2Tri;
209: dofA = dofAP2Tri;
210: }
211: break;
212: case 9: /* Quad */
213: if (order == 1) {
214: dofU = dofUP1Quad;
215: dofA = dofAP1Quad;
216: } else {
217: dofU = dofUP2Quad;
218: dofA = dofAP2Quad;
219: }
220: break;
221: case 15: /* Tet */
222: if (order == 1) {
223: dofU = dofUP1Tet;
224: dofA = dofAP1Tet;
225: } else {
226: dofU = dofUP2Tet;
227: dofA = dofAP2Tet;
228: }
229: break;
230: case 27: /* Hex */
231: if (order == 1) {
232: dofU = dofUP1Hex;
233: dofA = dofAP1Hex;
234: } else {
235: dofU = dofUP2Hex;
236: dofA = dofAP2Hex;
237: }
238: break;
239: default:
240: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Unknown element with closure size %" PetscInt_FMT, closureSize);
241: }
242: PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA));
244: for (cell = 0; cell < numCells; cell++) {
245: PetscInt *closure = NULL;
247: PetscCall(DMPlexGetTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure));
248: for (p = 0; p < closureSize; ++p) {
249: /* Find depth of p */
250: for (d = 0; d <= sdim; ++d) {
251: if ((closure[2 * p] >= pStartDepth[d]) && (closure[2 * p] < pEndDepth[d])) {
252: PetscCall(PetscSectionSetDof(section, closure[2 * p], dofU[d] + dofA[d] + dofS[d]));
253: PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldU, dofU[d]));
254: PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldA, dofA[d]));
255: PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldS, dofS[d]));
256: }
257: }
258: }
259: PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure));
260: }
261: PetscCall(ISRestoreIndices(cellIS, &cellID));
262: PetscCall(ISDestroy(&cellIS));
263: }
264: }
265: if (csIS) PetscCall(ISRestoreIndices(csIS, &csID));
266: PetscCall(ISDestroy(&csIS));
267: PetscCall(PetscSectionSetUp(section));
268: PetscCall(DMSetLocalSection(dm, section));
269: PetscCall(PetscObjectViewFromOptions((PetscObject)section, NULL, "-dm_section_view"));
270: PetscCall(PetscSectionDestroy(§ion));
272: {
273: PetscSF migrationSF;
274: PetscInt ovlp = 0;
275: PetscPartitioner part;
277: PetscCall(DMSetUseNatural(dm, PETSC_TRUE));
278: PetscCall(DMPlexGetPartitioner(dm, &part));
279: PetscCall(PetscPartitionerSetFromOptions(part));
280: PetscCall(DMPlexDistribute(dm, ovlp, &migrationSF, &pdm));
281: if (!pdm) pdm = dm;
282: /* Set the migrationSF is mandatory since useNatural = PETSC_TRUE */
283: if (migrationSF) {
284: PetscCall(DMPlexSetMigrationSF(pdm, migrationSF));
285: PetscCall(PetscSFDestroy(&migrationSF));
286: }
287: PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
288: }
290: /* Get DM and IS for each field of dm */
291: PetscCall(DMCreateSubDM(pdm, 1, &fieldU, &isU, &dmU));
292: PetscCall(DMCreateSubDM(pdm, 1, &fieldA, &isA, &dmA));
293: PetscCall(DMCreateSubDM(pdm, 1, &fieldS, &isS, &dmS));
294: PetscCall(DMCreateSubDM(pdm, 2, fieldUA, &isUA, &dmUA));
296: PetscCall(PetscMalloc1(2, &dmList));
297: dmList[0] = dmU;
298: dmList[1] = dmA;
300: PetscCall(DMCreateSuperDM(dmList, 2, NULL, &dmUA2));
301: PetscCall(PetscFree(dmList));
303: PetscCall(DMGetGlobalVector(pdm, &X));
304: PetscCall(DMGetGlobalVector(dmU, &U));
305: PetscCall(DMGetGlobalVector(dmA, &A));
306: PetscCall(DMGetGlobalVector(dmS, &S));
307: PetscCall(DMGetGlobalVector(dmUA, &UA));
308: PetscCall(DMGetGlobalVector(dmUA2, &UA2));
310: PetscCall(PetscObjectSetName((PetscObject)U, "U"));
311: PetscCall(PetscObjectSetName((PetscObject)A, "Alpha"));
312: PetscCall(PetscObjectSetName((PetscObject)S, "Sigma"));
313: PetscCall(PetscObjectSetName((PetscObject)UA, "UAlpha"));
314: PetscCall(PetscObjectSetName((PetscObject)UA2, "UAlpha2"));
315: PetscCall(VecSet(X, -111.));
317: /* Setting u to [x,y,z] and alpha to x^2+y^2+z^2 by writing in UAlpha then restricting to U and Alpha */
318: {
319: PetscSection sectionUA;
320: Vec UALoc;
321: PetscSection coordSection;
322: Vec coord;
323: PetscScalar *cval, *xyz;
324: PetscInt clSize, i, j;
326: PetscCall(DMGetLocalSection(dmUA, §ionUA));
327: PetscCall(DMGetLocalVector(dmUA, &UALoc));
328: PetscCall(VecGetArray(UALoc, &cval));
329: PetscCall(DMGetCoordinateSection(dmUA, &coordSection));
330: PetscCall(DMGetCoordinatesLocal(dmUA, &coord));
331: PetscCall(DMPlexGetChart(dmUA, &pStart, &pEnd));
333: for (p = pStart; p < pEnd; ++p) {
334: PetscInt dofUA, offUA;
336: PetscCall(PetscSectionGetDof(sectionUA, p, &dofUA));
337: if (dofUA > 0) {
338: xyz = NULL;
339: PetscCall(PetscSectionGetOffset(sectionUA, p, &offUA));
340: PetscCall(DMPlexVecGetClosure(dmUA, coordSection, coord, p, &clSize, &xyz));
341: cval[offUA + sdim] = 0.;
342: for (i = 0; i < sdim; ++i) {
343: cval[offUA + i] = 0;
344: for (j = 0; j < clSize / sdim; ++j) cval[offUA + i] += xyz[j * sdim + i];
345: cval[offUA + i] = cval[offUA + i] * sdim / clSize;
346: cval[offUA + sdim] += PetscSqr(cval[offUA + i]);
347: }
348: PetscCall(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, &clSize, &xyz));
349: }
350: }
351: PetscCall(VecRestoreArray(UALoc, &cval));
352: PetscCall(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA));
353: PetscCall(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA));
354: PetscCall(DMRestoreLocalVector(dmUA, &UALoc));
356: /* Update X */
357: PetscCall(VecISCopy(X, isUA, SCATTER_FORWARD, UA));
358: PetscCall(VecViewFromOptions(UA, NULL, "-ua_vec_view"));
360: /* Restrict to U and Alpha */
361: PetscCall(VecISCopy(X, isU, SCATTER_REVERSE, U));
362: PetscCall(VecISCopy(X, isA, SCATTER_REVERSE, A));
364: /* restrict to UA2 */
365: PetscCall(VecISCopy(X, isUA, SCATTER_REVERSE, UA2));
366: PetscCall(VecViewFromOptions(UA2, NULL, "-ua2_vec_view"));
367: }
369: {
370: Vec tmpVec;
371: PetscSection coordSection;
372: Vec coord;
373: PetscReal norm;
374: PetscReal time = 1.234;
376: /* Writing nodal variables to ExodusII file */
377: PetscCall(DMSetOutputSequenceNumber(dmU, 0, time));
378: PetscCall(DMSetOutputSequenceNumber(dmA, 0, time));
380: PetscCall(VecView(U, viewer));
381: PetscCall(VecView(A, viewer));
383: /* Saving U and Alpha in one shot.
384: For this, we need to cheat and change the Vec's name
385: Note that in the end we write variables one component at a time,
386: so that there is no real values in doing this
387: */
389: PetscCall(DMSetOutputSequenceNumber(dmUA, 1, time));
390: PetscCall(DMGetGlobalVector(dmUA, &tmpVec));
391: PetscCall(VecCopy(UA, tmpVec));
392: PetscCall(PetscObjectSetName((PetscObject)tmpVec, "U"));
393: PetscCall(VecView(tmpVec, viewer));
394: /* Reading nodal variables in Exodus file */
395: PetscCall(VecSet(tmpVec, -1000.0));
396: PetscCall(VecLoad(tmpVec, viewer));
397: PetscCall(VecAXPY(UA, -1.0, tmpVec));
398: PetscCall(VecNorm(UA, NORM_INFINITY, &norm));
399: PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UAlpha ||Vin - Vout|| = %g", (double)norm);
400: PetscCall(DMRestoreGlobalVector(dmUA, &tmpVec));
402: /* same thing with the UA2 Vec obtained from the superDM */
403: PetscCall(DMGetGlobalVector(dmUA2, &tmpVec));
404: PetscCall(VecCopy(UA2, tmpVec));
405: PetscCall(PetscObjectSetName((PetscObject)tmpVec, "U"));
406: PetscCall(DMSetOutputSequenceNumber(dmUA2, 2, time));
407: PetscCall(VecView(tmpVec, viewer));
408: /* Reading nodal variables in Exodus file */
409: PetscCall(VecSet(tmpVec, -1000.0));
410: PetscCall(VecLoad(tmpVec, viewer));
411: PetscCall(VecAXPY(UA2, -1.0, tmpVec));
412: PetscCall(VecNorm(UA2, NORM_INFINITY, &norm));
413: PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UAlpha2 ||Vin - Vout|| = %g", (double)norm);
414: PetscCall(DMRestoreGlobalVector(dmUA2, &tmpVec));
416: /* Building and saving Sigma
417: We set sigma_0 = rank (to see partitioning)
418: sigma_1 = cell set ID
419: sigma_2 = x_coordinate of the cell center of mass
420: */
421: PetscCall(DMGetCoordinateSection(dmS, &coordSection));
422: PetscCall(DMGetCoordinatesLocal(dmS, &coord));
423: PetscCall(DMGetLabelIdIS(dmS, "Cell Sets", &csIS));
424: PetscCall(DMGetLabelSize(dmS, "Cell Sets", &numCS));
425: PetscCall(ISGetIndices(csIS, &csID));
426: for (set = 0; set < numCS; ++set) {
427: /* We know that all cells in a cell set have the same type, so we can dimension cval and xyz once for each cell set */
428: IS cellIS;
429: const PetscInt *cellID;
430: PetscInt numCells, cell;
431: PetscScalar *cval = NULL, *xyz = NULL;
432: PetscInt clSize, cdimCoord, c;
434: PetscCall(DMGetStratumIS(dmS, "Cell Sets", csID[set], &cellIS));
435: PetscCall(ISGetIndices(cellIS, &cellID));
436: PetscCall(ISGetSize(cellIS, &numCells));
437: for (cell = 0; cell < numCells; cell++) {
438: PetscCall(DMPlexVecGetClosure(dmS, NULL, S, cellID[cell], &clSize, &cval));
439: PetscCall(DMPlexVecGetClosure(dmS, coordSection, coord, cellID[cell], &cdimCoord, &xyz));
440: cval[0] = rank;
441: cval[1] = csID[set];
442: cval[2] = 0.;
443: for (c = 0; c < cdimCoord / sdim; c++) cval[2] += xyz[c * sdim];
444: cval[2] = cval[2] * sdim / cdimCoord;
445: PetscCall(DMPlexVecSetClosure(dmS, NULL, S, cellID[cell], cval, INSERT_ALL_VALUES));
446: }
447: PetscCall(DMPlexVecRestoreClosure(dmS, NULL, S, cellID[0], &clSize, &cval));
448: PetscCall(DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID[0], NULL, &xyz));
449: PetscCall(ISRestoreIndices(cellIS, &cellID));
450: PetscCall(ISDestroy(&cellIS));
451: }
452: PetscCall(ISRestoreIndices(csIS, &csID));
453: PetscCall(ISDestroy(&csIS));
454: PetscCall(VecViewFromOptions(S, NULL, "-s_vec_view"));
456: /* Writing zonal variables in Exodus file */
457: PetscCall(DMSetOutputSequenceNumber(dmS, 0, time));
458: PetscCall(VecView(S, viewer));
460: /* Reading zonal variables in Exodus file */
461: PetscCall(DMGetGlobalVector(dmS, &tmpVec));
462: PetscCall(VecSet(tmpVec, -1000.0));
463: PetscCall(PetscObjectSetName((PetscObject)tmpVec, "Sigma"));
464: PetscCall(VecLoad(tmpVec, viewer));
465: PetscCall(VecAXPY(S, -1.0, tmpVec));
466: PetscCall(VecNorm(S, NORM_INFINITY, &norm));
467: PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Sigma ||Vin - Vout|| = %g", (double)norm);
468: PetscCall(DMRestoreGlobalVector(dmS, &tmpVec));
469: }
470: PetscCall(PetscViewerDestroy(&viewer));
472: PetscCall(DMRestoreGlobalVector(dmUA2, &UA2));
473: PetscCall(DMRestoreGlobalVector(dmUA, &UA));
474: PetscCall(DMRestoreGlobalVector(dmS, &S));
475: PetscCall(DMRestoreGlobalVector(dmA, &A));
476: PetscCall(DMRestoreGlobalVector(dmU, &U));
477: PetscCall(DMRestoreGlobalVector(pdm, &X));
478: PetscCall(DMDestroy(&dmU));
479: PetscCall(ISDestroy(&isU));
480: PetscCall(DMDestroy(&dmA));
481: PetscCall(ISDestroy(&isA));
482: PetscCall(DMDestroy(&dmS));
483: PetscCall(ISDestroy(&isS));
484: PetscCall(DMDestroy(&dmUA));
485: PetscCall(ISDestroy(&isUA));
486: PetscCall(DMDestroy(&dmUA2));
487: PetscCall(DMDestroy(&pdm));
488: if (size > 1) PetscCall(DMDestroy(&dm));
489: PetscCall(PetscFree2(pStartDepth, pEndDepth));
490: PetscCall(PetscFinalize());
491: return 0;
492: }
494: /*TEST
496: build:
497: requires: exodusii pnetcdf !complex
498: # 2D seq
499: test:
500: suffix: 0
501: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 1
502: #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
503: test:
504: suffix: 1
505: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 1
506: test:
507: suffix: 2
508: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 1
509: #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
510: test:
511: suffix: 3
512: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 2
513: test:
514: suffix: 4
515: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 2
516: test:
517: suffix: 5
518: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 2
520: # 2D par
521: test:
522: suffix: 6
523: nsize: 2
524: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 1
525: #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
526: test:
527: suffix: 7
528: nsize: 2
529: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 1
530: test:
531: suffix: 8
532: nsize: 2
533: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 1
534: #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name
535: test:
536: suffix: 9
537: nsize: 2
538: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 2
539: test:
540: suffix: 10
541: nsize: 2
542: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 2
543: test:
544: # Something is now broken with parallel read/write for whatever shape H is
545: TODO: broken
546: suffix: 11
547: nsize: 2
548: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 2
550: #3d seq
551: test:
552: suffix: 12
553: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 1
554: test:
555: suffix: 13
556: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 1
557: #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
558: test:
559: suffix: 14
560: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 2
561: test:
562: suffix: 15
563: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 2
564: #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
565: #3d par
566: test:
567: suffix: 16
568: nsize: 2
569: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 1
570: test:
571: suffix: 17
572: nsize: 2
573: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 1
574: #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
575: test:
576: suffix: 18
577: nsize: 2
578: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 2
579: test:
580: suffix: 19
581: nsize: 2
582: args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 2
583: #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
585: TEST*/