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. 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 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 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 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 ${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*/ 143