Actual source code: sectionhdf5.c

      1: #include <petsc/private/sectionimpl.h>
      2: #include <petscsf.h>
      3: #include <petscis.h>
      4: #include <petscviewerhdf5.h>
      5: #include <petsclayouthdf5.h>

      7: static PetscErrorCode PetscSectionView_HDF5_SingleField(PetscSection s, PetscViewer viewer)
      8: {
      9:   MPI_Comm  comm;
     10:   PetscInt  pStart, pEnd, p, n;
     11:   PetscBool hasConstraints, includesConstraints;
     12:   IS        dofIS, offIS, cdofIS, coffIS, cindIS;
     13:   PetscInt *dofs, *offs, *cdofs, *coffs, *cinds, dof, cdof, m, moff, i;

     15:   PetscFunctionBegin;
     16:   PetscCall(PetscObjectGetComm((PetscObject)s, &comm));
     17:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
     18:   hasConstraints = (s->bc) ? PETSC_TRUE : PETSC_FALSE;
     19:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hasConstraints, 1, MPIU_BOOL, MPI_LOR, comm));
     20:   for (p = pStart, n = 0, m = 0; p < pEnd; ++p) {
     21:     PetscCall(PetscSectionGetDof(s, p, &dof));
     22:     if (dof >= 0) {
     23:       if (hasConstraints) {
     24:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
     25:         m += cdof;
     26:       }
     27:       n++;
     28:     }
     29:   }
     30:   PetscCall(PetscMalloc1(n, &dofs));
     31:   PetscCall(PetscMalloc1(n, &offs));
     32:   if (hasConstraints) {
     33:     PetscCall(PetscMalloc1(n, &cdofs));
     34:     PetscCall(PetscMalloc1(n, &coffs));
     35:     PetscCall(PetscMalloc1(m, &cinds));
     36:   }
     37:   for (p = pStart, n = 0, m = 0; p < pEnd; ++p) {
     38:     PetscCall(PetscSectionGetDof(s, p, &dof));
     39:     if (dof >= 0) {
     40:       dofs[n] = dof;
     41:       PetscCall(PetscSectionGetOffset(s, p, &offs[n]));
     42:       if (hasConstraints) {
     43:         const PetscInt *cpinds;

     45:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
     46:         PetscCall(PetscSectionGetConstraintIndices(s, p, &cpinds));
     47:         cdofs[n] = cdof;
     48:         coffs[n] = m;
     49:         for (i = 0; i < cdof; ++i) cinds[m++] = cpinds[i];
     50:       }
     51:       n++;
     52:     }
     53:   }
     54:   if (hasConstraints) {
     55:     PetscCallMPI(MPI_Scan(&m, &moff, 1, MPIU_INT, MPI_SUM, comm));
     56:     moff -= m;
     57:     for (p = 0; p < n; ++p) coffs[p] += moff;
     58:   }
     59:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "hasConstraints", PETSC_BOOL, (void *)&hasConstraints));
     60:   PetscCall(PetscSectionGetIncludesConstraints(s, &includesConstraints));
     61:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "includesConstraints", PETSC_BOOL, (void *)&includesConstraints));
     62:   PetscCall(ISCreateGeneral(comm, n, dofs, PETSC_OWN_POINTER, &dofIS));
     63:   PetscCall(PetscObjectSetName((PetscObject)dofIS, "atlasDof"));
     64:   PetscCall(ISView(dofIS, viewer));
     65:   PetscCall(ISDestroy(&dofIS));
     66:   PetscCall(ISCreateGeneral(comm, n, offs, PETSC_OWN_POINTER, &offIS));
     67:   PetscCall(PetscObjectSetName((PetscObject)offIS, "atlasOff"));
     68:   PetscCall(ISView(offIS, viewer));
     69:   PetscCall(ISDestroy(&offIS));
     70:   if (hasConstraints) {
     71:     PetscCall(PetscViewerHDF5PushGroup(viewer, "bc"));
     72:     PetscCall(ISCreateGeneral(comm, n, cdofs, PETSC_OWN_POINTER, &cdofIS));
     73:     PetscCall(PetscObjectSetName((PetscObject)cdofIS, "atlasDof"));
     74:     PetscCall(ISView(cdofIS, viewer));
     75:     PetscCall(ISDestroy(&cdofIS));
     76:     PetscCall(ISCreateGeneral(comm, n, coffs, PETSC_OWN_POINTER, &coffIS));
     77:     PetscCall(PetscObjectSetName((PetscObject)coffIS, "atlasOff"));
     78:     PetscCall(ISView(coffIS, viewer));
     79:     PetscCall(ISDestroy(&coffIS));
     80:     PetscCall(PetscViewerHDF5PopGroup(viewer));
     81:     PetscCall(ISCreateGeneral(comm, m, cinds, PETSC_OWN_POINTER, &cindIS));
     82:     PetscCall(PetscObjectSetName((PetscObject)cindIS, "bcIndices"));
     83:     PetscCall(ISView(cindIS, viewer));
     84:     PetscCall(ISDestroy(&cindIS));
     85:   }
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: PetscErrorCode PetscSectionView_HDF5_Internal(PetscSection s, PetscViewer viewer)
     90: {
     91:   PetscInt numFields, f;

     93:   PetscFunctionBegin;
     94:   PetscCall(PetscViewerHDF5PushGroup(viewer, "section"));
     95:   PetscCall(PetscSectionGetNumFields(s, &numFields));
     96:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "numFields", PETSC_INT, (void *)&numFields));
     97:   PetscCall(PetscSectionView_HDF5_SingleField(s, viewer));
     98:   for (f = 0; f < numFields; ++f) {
     99:     char        fname[PETSC_MAX_PATH_LEN];
    100:     const char *fieldName;
    101:     PetscInt    fieldComponents, c;

    103:     PetscCall(PetscSNPrintf(fname, sizeof(fname), "field%" PetscInt_FMT, f));
    104:     PetscCall(PetscViewerHDF5PushGroup(viewer, fname));
    105:     PetscCall(PetscSectionGetFieldName(s, f, &fieldName));
    106:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "fieldName", PETSC_STRING, fieldName));
    107:     PetscCall(PetscSectionGetFieldComponents(s, f, &fieldComponents));
    108:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "fieldComponents", PETSC_INT, (void *)&fieldComponents));
    109:     for (c = 0; c < fieldComponents; ++c) {
    110:       char        cname[PETSC_MAX_PATH_LEN];
    111:       const char *componentName;

    113:       PetscCall(PetscSNPrintf(cname, sizeof(cname), "component%" PetscInt_FMT, c));
    114:       PetscCall(PetscViewerHDF5PushGroup(viewer, cname));
    115:       PetscCall(PetscSectionGetComponentName(s, f, c, &componentName));
    116:       PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "componentName", PETSC_STRING, componentName));
    117:       PetscCall(PetscViewerHDF5PopGroup(viewer));
    118:     }
    119:     PetscCall(PetscSectionView_HDF5_SingleField(s->field[f], viewer));
    120:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    121:   }
    122:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    123:   PetscFunctionReturn(PETSC_SUCCESS);
    124: }

    126: static PetscErrorCode PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(PetscSection s, IS cindIS, IS coffIS)
    127: {
    128:   MPI_Comm        comm;
    129:   PetscInt        pStart, pEnd, p, M, m, i, cdof;
    130:   const PetscInt *data;
    131:   PetscInt       *cinds;
    132:   const PetscInt *coffs;
    133:   PetscInt       *coffsets;
    134:   PetscSF         sf;
    135:   PetscLayout     layout;

    137:   PetscFunctionBegin;
    138:   PetscCall(PetscObjectGetComm((PetscObject)s, &comm));
    139:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
    140:   PetscCall(ISGetSize(cindIS, &M));
    141:   PetscCall(ISGetLocalSize(cindIS, &m));
    142:   PetscCall(PetscMalloc1(m, &coffsets));
    143:   PetscCall(ISGetIndices(coffIS, &coffs));
    144:   for (p = pStart, m = 0; p < pEnd; ++p) {
    145:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
    146:     for (i = 0; i < cdof; ++i) coffsets[m++] = coffs[p - pStart] + i;
    147:   }
    148:   PetscCall(ISRestoreIndices(coffIS, &coffs));
    149:   PetscCall(PetscSFCreate(comm, &sf));
    150:   PetscCall(PetscLayoutCreate(comm, &layout));
    151:   PetscCall(PetscLayoutSetSize(layout, M));
    152:   PetscCall(PetscLayoutSetLocalSize(layout, m));
    153:   PetscCall(PetscLayoutSetBlockSize(layout, 1));
    154:   PetscCall(PetscLayoutSetUp(layout));
    155:   PetscCall(PetscSFSetGraphLayout(sf, layout, m, NULL, PETSC_OWN_POINTER, coffsets));
    156:   PetscCall(PetscLayoutDestroy(&layout));
    157:   PetscCall(PetscFree(coffsets));
    158:   PetscCall(PetscMalloc1(m, &cinds));
    159:   PetscCall(ISGetIndices(cindIS, &data));
    160:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, data, cinds, MPI_REPLACE));
    161:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, data, cinds, MPI_REPLACE));
    162:   PetscCall(ISRestoreIndices(cindIS, &data));
    163:   PetscCall(PetscSFDestroy(&sf));
    164:   PetscCall(PetscSectionSetUpBC(s));
    165:   for (p = pStart, m = 0; p < pEnd; ++p) {
    166:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
    167:     PetscCall(PetscSectionSetConstraintIndices(s, p, PetscSafePointerPlusOffset(cinds, m)));
    168:     m += cdof;
    169:   }
    170:   PetscCall(PetscFree(cinds));
    171:   PetscFunctionReturn(PETSC_SUCCESS);
    172: }

    174: static PetscErrorCode PetscSectionLoad_HDF5_SingleField(PetscSection s, PetscViewer viewer)
    175: {
    176:   MPI_Comm comm;
    177:   PetscInt pStart, pEnd, p, N, n, M, m;
    178: #if defined(PETSC_USE_DEBUG)
    179:   PetscInt N1, M1;
    180: #endif
    181:   PetscBool       hasConstraints, includesConstraints;
    182:   IS              dofIS, offIS, cdofIS, coffIS, cindIS;
    183:   const PetscInt *dofs, *offs, *cdofs;
    184:   PetscLayout     map;

    186:   PetscFunctionBegin;
    187:   PetscCall(PetscObjectGetComm((PetscObject)s, &comm));
    188:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "includesConstraints", PETSC_BOOL, NULL, (void *)&includesConstraints));
    189:   PetscCall(PetscSectionSetIncludesConstraints(s, includesConstraints));
    190:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
    191:   n = pEnd - pStart;
    192: #if defined(PETSC_USE_DEBUG)
    193:   PetscCallMPI(MPIU_Allreduce(&n, &N1, 1, MPIU_INT, MPI_SUM, comm));
    194: #endif
    195:   PetscCall(ISCreate(comm, &dofIS));
    196:   PetscCall(PetscObjectSetName((PetscObject)dofIS, "atlasDof"));
    197:   PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N));
    198: #if defined(PETSC_USE_DEBUG)
    199:   PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->atlasDof: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, N1, N, n);
    200: #endif
    201:   PetscCall(ISGetLayout(dofIS, &map));
    202:   PetscCall(PetscLayoutSetSize(map, N));
    203:   PetscCall(PetscLayoutSetLocalSize(map, n));
    204:   PetscCall(ISLoad(dofIS, viewer));
    205:   PetscCall(ISCreate(comm, &offIS));
    206:   PetscCall(PetscObjectSetName((PetscObject)offIS, "atlasOff"));
    207:   PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasOff", NULL, &N));
    208: #if defined(PETSC_USE_DEBUG)
    209:   PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->atlasOff: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, N1, N, n);
    210: #endif
    211:   PetscCall(ISGetLayout(offIS, &map));
    212:   PetscCall(PetscLayoutSetSize(map, N));
    213:   PetscCall(PetscLayoutSetLocalSize(map, n));
    214:   PetscCall(ISLoad(offIS, viewer));
    215:   PetscCall(ISGetIndices(dofIS, &dofs));
    216:   PetscCall(ISGetIndices(offIS, &offs));
    217:   for (p = pStart, n = 0; p < pEnd; ++p, ++n) {
    218:     PetscCall(PetscSectionSetDof(s, p, dofs[n]));
    219:     PetscCall(PetscSectionSetOffset(s, p, offs[n]));
    220:   }
    221:   PetscCall(ISRestoreIndices(dofIS, &dofs));
    222:   PetscCall(ISRestoreIndices(offIS, &offs));
    223:   PetscCall(ISDestroy(&dofIS));
    224:   PetscCall(ISDestroy(&offIS));
    225:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "hasConstraints", PETSC_BOOL, NULL, (void *)&hasConstraints));
    226:   if (hasConstraints) {
    227:     PetscCall(PetscViewerHDF5PushGroup(viewer, "bc"));
    228:     PetscCall(ISCreate(comm, &cdofIS));
    229:     PetscCall(PetscObjectSetName((PetscObject)cdofIS, "atlasDof"));
    230:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N));
    231: #if defined(PETSC_USE_DEBUG)
    232:     PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->bc->atlasDof: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, N1, N, n);
    233: #endif
    234:     PetscCall(ISGetLayout(cdofIS, &map));
    235:     PetscCall(PetscLayoutSetSize(map, N));
    236:     PetscCall(PetscLayoutSetLocalSize(map, n));
    237:     PetscCall(ISLoad(cdofIS, viewer));
    238:     PetscCall(ISGetIndices(cdofIS, &cdofs));
    239:     for (p = pStart, n = 0; p < pEnd; ++p, ++n) PetscCall(PetscSectionSetConstraintDof(s, p, cdofs[n]));
    240:     PetscCall(ISRestoreIndices(cdofIS, &cdofs));
    241:     PetscCall(ISDestroy(&cdofIS));
    242:     PetscCall(ISCreate(comm, &coffIS));
    243:     PetscCall(PetscObjectSetName((PetscObject)coffIS, "atlasOff"));
    244:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasOff", NULL, &N));
    245: #if defined(PETSC_USE_DEBUG)
    246:     PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->bc->atlasOff: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, N1, N, n);
    247: #endif
    248:     PetscCall(ISGetLayout(coffIS, &map));
    249:     PetscCall(PetscLayoutSetSize(map, N));
    250:     PetscCall(PetscLayoutSetLocalSize(map, n));
    251:     PetscCall(ISLoad(coffIS, viewer));
    252:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    253:     PetscCall(ISCreate(comm, &cindIS));
    254:     PetscCall(PetscObjectSetName((PetscObject)cindIS, "bcIndices"));
    255:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "bcIndices", NULL, &M));
    256:     if (!s->bc) m = 0;
    257:     else PetscCall(PetscSectionGetStorageSize(s->bc, &m));
    258: #if defined(PETSC_USE_DEBUG)
    259:     PetscCallMPI(MPIU_Allreduce(&m, &M1, 1, MPIU_INT, MPI_SUM, comm));
    260:     PetscCheck(M1 == M, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->bcIndices: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, M1, M, m);
    261: #endif
    262:     PetscCall(ISGetLayout(cindIS, &map));
    263:     PetscCall(PetscLayoutSetSize(map, M));
    264:     PetscCall(PetscLayoutSetLocalSize(map, m));
    265:     PetscCall(ISLoad(cindIS, viewer));
    266:     PetscCall(PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(s, cindIS, coffIS));
    267:     PetscCall(ISDestroy(&coffIS));
    268:     PetscCall(ISDestroy(&cindIS));
    269:   }
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: PetscErrorCode PetscSectionLoad_HDF5_Internal(PetscSection s, PetscViewer viewer)
    274: {
    275:   MPI_Comm comm;
    276:   PetscInt N, n, numFields, f;

    278:   PetscFunctionBegin;
    279:   PetscCall(PetscObjectGetComm((PetscObject)s, &comm));
    280:   PetscCall(PetscViewerHDF5PushGroup(viewer, "section"));
    281:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "numFields", PETSC_INT, NULL, (void *)&numFields));
    282:   if (s->pStart < 0 && s->pEnd < 0) n = PETSC_DECIDE;
    283:   else {
    284:     PetscCheck(s->pStart == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "s->pStart must be 0 (got %" PetscInt_FMT ")", s->pStart);
    285:     PetscCheck(s->pEnd >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "s->pEnd must be >= 0, (got %" PetscInt_FMT ")", s->pEnd);
    286:     n = s->pEnd;
    287:   }
    288:   if (numFields > 0) PetscCall(PetscSectionSetNumFields(s, numFields));
    289:   PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N));
    290:   if (n == PETSC_DECIDE) PetscCall(PetscSplitOwnership(comm, &n, &N));
    291:   PetscCall(PetscSectionSetChart(s, 0, n));
    292:   PetscCall(PetscSectionLoad_HDF5_SingleField(s, viewer));
    293:   for (f = 0; f < numFields; ++f) {
    294:     char     fname[PETSC_MAX_PATH_LEN];
    295:     char    *fieldName;
    296:     PetscInt fieldComponents, c;

    298:     PetscCall(PetscSNPrintf(fname, sizeof(fname), "field%" PetscInt_FMT, f));
    299:     PetscCall(PetscViewerHDF5PushGroup(viewer, fname));
    300:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldName", PETSC_STRING, NULL, &fieldName));
    301:     PetscCall(PetscSectionSetFieldName(s, f, fieldName));
    302:     PetscCall(PetscFree(fieldName));
    303:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldComponents", PETSC_INT, NULL, (void *)&fieldComponents));
    304:     PetscCall(PetscSectionSetFieldComponents(s, f, fieldComponents));
    305:     for (c = 0; c < fieldComponents; ++c) {
    306:       char  cname[PETSC_MAX_PATH_LEN];
    307:       char *componentName;

    309:       PetscCall(PetscSNPrintf(cname, sizeof(cname), "component%" PetscInt_FMT, c));
    310:       PetscCall(PetscViewerHDF5PushGroup(viewer, cname));
    311:       PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "componentName", PETSC_STRING, NULL, &componentName));
    312:       PetscCall(PetscSectionSetComponentName(s, f, c, componentName));
    313:       PetscCall(PetscFree(componentName));
    314:       PetscCall(PetscViewerHDF5PopGroup(viewer));
    315:     }
    316:     PetscCall(PetscSectionLoad_HDF5_SingleField(s->field[f], viewer));
    317:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    318:   }
    319:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    320:   PetscFunctionReturn(PETSC_SUCCESS);
    321: }