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 initital
  6: // CGNS file to support performance benchmarking.

  8: #include <petscdmplex.h>
  9: #include <petscviewerhdf5.h>
 10: #define EX "ex16.c"

 12: typedef struct {
 13:   char infile[PETSC_MAX_PATH_LEN]; /* Input mesh filename */
 14: } AppCtx;

 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: }

 27: // @brief Create DM from CGNS file and setup PetscFE to VecLoad solution from that file
 28: PetscErrorCode ReadCGNSDM(MPI_Comm comm, const char filename[], DM *dm)
 29: {
 30:   PetscInt degree;

 32:   PetscFunctionBeginUser;
 33:   PetscCall(DMPlexCreateFromFile(comm, filename, "ex16_plex", PETSC_TRUE, dm));
 34:   PetscCall(DMSetFromOptions(*dm));
 35:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));

 37:   { // Get degree of the natural section
 38:     PetscFE        fe_natural;
 39:     PetscDualSpace dual_space_natural;

 41:     PetscCall(DMGetField(*dm, 0, NULL, (PetscObject *)&fe_natural));
 42:     PetscCall(PetscFEGetDualSpace(fe_natural, &dual_space_natural));
 43:     PetscCall(PetscDualSpaceGetOrder(dual_space_natural, &degree));
 44:     PetscCall(DMClearFields(*dm));
 45:     PetscCall(DMSetLocalSection(*dm, NULL));
 46:   }

 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;

 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:   }

 74:   // Set section component names, used when writing out CGNS files
 75:   PetscSection section;
 76:   PetscCall(DMGetLocalSection(*dm, &section));
 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: }

 86: int main(int argc, char **argv)
 87: {
 88:   AppCtx      user;
 89:   MPI_Comm    comm;
 90:   const char *infilename;

 92:   PetscFunctionBeginUser;
 93:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 94:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
 95:   infilename = user.infile;

 97:   DM          dm;
 98:   Vec         V;
 99:   PetscViewer viewer;
100:   const char *name;
101:   PetscReal   time;
102:   PetscBool   set;
103:   comm = PETSC_COMM_WORLD;

105:   // Load DM from CGNS file
106:   PetscCall(ReadCGNSDM(comm, infilename, &dm));
107:   PetscCall(DMSetOptionsPrefix(dm, "loaded_"));
108:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

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));

120:   // Write loaded solution (e.g. example in TEST below is to CGNS file)
121:   PetscCall(VecViewFromOptions(V, NULL, "-vec_view"));

123:   PetscCall(DMRestoreGlobalVector(dm, &V));
124:   PetscCall(DMDestroy(&dm));

126:   PetscCall(PetscFinalize());
127:   return 0;
128: }

130: /*TEST
131:   build:
132:     requires: cgns
133:   testset:
134:     suffix: cgns
135:     requires: !complex
136:     nsize: 4
137:     args: -infile ${wPETSC_DIR}/share/petsc/datafiles/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*/