1 #include "petscsf.h"
2 static char help[] = "Simple demonstration of CGNS parallel load-save including data\n\n";
3 // As this is a tutorial that is intended to be an easy starting point feel free to make new
4 // example files that extend this but please keep this one simple.
5 // In subsequent examples we will also provide tools to generate an arbitrary size initial
6 // CGNS file to support performance benchmarking.
7
8 #include <petscdmplex.h>
9 #include <petscviewerhdf5.h>
10 #define EX "ex16.c"
11
12 typedef struct {
13 char infile[PETSC_MAX_PATH_LEN]; /* Input mesh filename */
14 } AppCtx;
15
ProcessOptions(MPI_Comm comm,AppCtx * options)16 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
17 {
18 PetscFunctionBeginUser;
19 options->infile[0] = '\0';
20 PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
21 PetscCall(PetscOptionsString("-infile", "The input CGNS file", EX, options->infile, options->infile, sizeof(options->infile), NULL));
22 PetscOptionsEnd();
23 PetscCheck(options->infile[0], comm, PETSC_ERR_USER_INPUT, "-infile needs to be specified");
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
27 // @brief Create DM from CGNS file and setup PetscFE to VecLoad solution from that file
ReadCGNSDM(MPI_Comm comm,const char filename[],DM * dm)28 PetscErrorCode ReadCGNSDM(MPI_Comm comm, const char filename[], DM *dm)
29 {
30 PetscInt degree;
31
32 PetscFunctionBeginUser;
33 PetscCall(DMPlexCreateFromFile(comm, filename, "ex16_plex", PETSC_TRUE, dm));
34 PetscCall(DMSetFromOptions(*dm));
35 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
36
37 { // Get degree of the natural section
38 PetscFE fe_natural;
39 PetscDualSpace dual_space_natural;
40
41 PetscCall(DMGetField(*dm, 0, NULL, (PetscObject *)&fe_natural));
42 PetscCall(PetscFEGetDualSpace(fe_natural, &dual_space_natural));
43 PetscCall(PetscDualSpaceGetOrder(dual_space_natural, °ree));
44 PetscCall(DMClearFields(*dm));
45 PetscCall(DMSetLocalSection(*dm, NULL));
46 }
47
48 { // Setup fe to load in the initial condition data
49 PetscFE fe;
50 PetscInt dim, cStart, cEnd;
51 PetscInt ctInt, mincti, maxcti;
52 DMPolytopeType dm_polytope, cti;
53
54 PetscCall(DMGetDimension(*dm, &dim));
55 // Limiting to single topology in this simple example
56 PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd));
57 PetscCall(DMPlexGetCellType(*dm, cStart, &dm_polytope));
58 for (PetscInt i = cStart + 1; i < cEnd; i++) {
59 PetscCall(DMPlexGetCellType(*dm, i, &cti));
60 PetscCheck(cti == dm_polytope, comm, PETSC_ERR_RETURN, "Multi-topology not yet supported in this example!");
61 }
62 ctInt = cti;
63 PetscCallMPI(MPIU_Allreduce(&ctInt, &maxcti, 1, MPIU_INT, MPI_MAX, comm));
64 PetscCallMPI(MPIU_Allreduce(&ctInt, &mincti, 1, MPIU_INT, MPI_MIN, comm));
65 PetscCheck(mincti == maxcti, comm, PETSC_ERR_RETURN, "Multi-topology not yet supported in this example!");
66 PetscCall(PetscPrintf(comm, "Mesh confirmed to be single topology degree %" PetscInt_FMT " %s\n", degree, DMPolytopeTypes[cti]));
67 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 5, dm_polytope, degree, PETSC_DETERMINE, &fe));
68 PetscCall(PetscObjectSetName((PetscObject)fe, "FE for VecLoad"));
69 PetscCall(DMAddField(*dm, NULL, (PetscObject)fe));
70 PetscCall(DMCreateDS(*dm));
71 PetscCall(PetscFEDestroy(&fe));
72 }
73
74 // Set section component names, used when writing out CGNS files
75 PetscSection section;
76 PetscCall(DMGetLocalSection(*dm, §ion));
77 PetscCall(PetscSectionSetFieldName(section, 0, ""));
78 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Pressure"));
79 PetscCall(PetscSectionSetComponentName(section, 0, 1, "VelocityX"));
80 PetscCall(PetscSectionSetComponentName(section, 0, 2, "VelocityY"));
81 PetscCall(PetscSectionSetComponentName(section, 0, 3, "VelocityZ"));
82 PetscCall(PetscSectionSetComponentName(section, 0, 4, "Temperature"));
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85
main(int argc,char ** argv)86 int main(int argc, char **argv)
87 {
88 AppCtx user;
89 MPI_Comm comm;
90 const char *infilename;
91
92 PetscFunctionBeginUser;
93 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
94 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
95 infilename = user.infile;
96
97 DM dm;
98 Vec V;
99 PetscViewer viewer;
100 const char *name;
101 PetscReal time;
102 PetscBool set;
103 comm = PETSC_COMM_WORLD;
104
105 // Load DM from CGNS file
106 PetscCall(ReadCGNSDM(comm, infilename, &dm));
107 PetscCall(DMSetOptionsPrefix(dm, "loaded_"));
108 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
109
110 // Load solution from CGNS file
111 PetscCall(PetscViewerCGNSOpen(comm, infilename, FILE_MODE_READ, &viewer));
112 PetscCall(DMGetGlobalVector(dm, &V));
113 PetscCall(PetscViewerCGNSSetSolutionIndex(viewer, 1));
114 PetscCall(PetscViewerCGNSGetSolutionName(viewer, &name));
115 PetscCall(PetscViewerCGNSGetSolutionTime(viewer, &time, &set));
116 PetscCall(PetscPrintf(comm, "Solution Name: %s, and time %g\n", name, (double)time));
117 PetscCall(VecLoad(V, viewer));
118 PetscCall(PetscViewerDestroy(&viewer));
119
120 // Write loaded solution (e.g. example in TEST below is to CGNS file)
121 PetscCall(VecViewFromOptions(V, NULL, "-vec_view"));
122
123 PetscCall(DMRestoreGlobalVector(dm, &V));
124 PetscCall(DMDestroy(&dm));
125
126 PetscCall(PetscFinalize());
127 return 0;
128 }
129
130 /*TEST
131 build:
132 requires: cgns
133 testset:
134 suffix: cgns
135 requires: !complex
136 nsize: 4
137 args: -infile ${DATAFILESPATH}/meshes/2x2x2_Q3_wave.cgns
138 args: -dm_plex_cgns_parallel -loaded_dm_view
139 test:
140 suffix: simple
141 args: -vec_view cgns:2x2x2_Q3Vecview.cgns
142 TEST*/
143