1: #include "petsc/finclude/petscdmplex.h"
  2: program ex26f90
  3:   use petscdmplex
  4:   implicit none
  5: #include "exodusII.inc"

  7:   type(tDM)                           :: dm, dmU, dmA, dmS, dmUA, dmUA2, pDM
  8:   type(tDM), dimension(:), pointer    :: dmList
  9:   type(tVec)                          :: X, U, A, S, UA, UA2
 10:   type(tIS)                           :: isU, isA, isS, isUA
 11:   type(tPetscSection)                 :: section
 12:   PetscInt, parameter                 :: fieldU = 0, fieldS = 1, fieldA = 2
 13:   PetscInt, parameter, dimension(2)   :: fieldUA = [0, 2]
 14:   character(len=PETSC_MAX_PATH_LEN)   :: ifilename, ofilename, IOBuffer
 15:   integer                             :: exoid = -1
 16:   type(tIS)                           :: csIS
 17:   PetscInt, dimension(:), pointer     :: csID
 18:   PetscInt, dimension(:), pointer     :: pStartDepth, pEndDepth
 19:   PetscInt                            :: order = 1
 20:   integer                             :: i
 21:   PetscInt                            :: sdim, d, pStart, pEnd, p, numCS, set, j, k
 22:   PetscMPIInt                         :: rank, numProc
 23:   PetscBool                           :: flg
 24:   PetscErrorCode                      :: ierr
 25:   MPIU_Comm                           :: comm
 26:   type(tPetscViewer)                  :: viewer

 28:   character(len=MXSTLN)               :: sJunk
 29:   PetscInt, parameter                 :: numstep = 3
 30:   PetscInt                            :: step
 31:   integer                             :: numNodalVar, numZonalVar
 32:   character(len=MXNAME), dimension(4) :: nodalVarName2D = ["U_x  ", &
 33:                                                            "U_y  ", &
 34:                                                            "Alpha", &
 35:                                                            "Beta "]
 36:   character(len=MXNAME), dimension(3) :: zonalVarName2D = ["Sigma_11", &
 37:                                                            "Sigma_12", &
 38:                                                            "Sigma_22"]
 39:   character(len=MXNAME), dimension(5) :: nodalVarName3D = ["U_x  ", &
 40:                                                            "U_y  ", &
 41:                                                            "U_z  ", &
 42:                                                            "Alpha", &
 43:                                                            "Beta "]
 44:   character(len=MXNAME), dimension(6) :: zonalVarName3D = ["Sigma_11", &
 45:                                                            "Sigma_22", &
 46:                                                            "Sigma_33", &
 47:                                                            "Sigma_23", &
 48:                                                            "Sigma_13", &
 49:                                                            "Sigma_12"]
 50:   logical, dimension(:, :), pointer   :: truthtable

 52:   type(tIS)                           :: cellIS
 53:   PetscInt, dimension(:), pointer     :: cellID
 54:   PetscInt                            :: numCells, cell, closureSize
 55:   PetscInt, dimension(:), pointer     :: closureA, closure

 57:   type(tPetscSection)                 :: sectionUA, coordSection
 58:   type(tVec)                          :: UALoc, coord
 59:   PetscScalar, dimension(:), pointer  :: cval, xyz
 60:   PetscInt                            :: dofUA, offUA, c

 62:   ! dof layout ordered by increasing height in the DAG: cell, face, edge, vertex
 63:   PetscInt, dimension(3), target      :: dofS2D = [0, 0, 3]
 64:   PetscInt, dimension(3), target      :: dofUP1Tri = [2, 0, 0]
 65:   PetscInt, dimension(3), target      :: dofAP1Tri = [1, 0, 0]
 66:   PetscInt, dimension(3), target      :: dofUP2Tri = [2, 2, 0]
 67:   PetscInt, dimension(3), target      :: dofAP2Tri = [1, 1, 0]
 68:   PetscInt, dimension(3), target      :: dofUP1Quad = [2, 0, 0]
 69:   PetscInt, dimension(3), target      :: dofAP1Quad = [1, 0, 0]
 70:   PetscInt, dimension(3), target      :: dofUP2Quad = [2, 2, 2]
 71:   PetscInt, dimension(3), target      :: dofAP2Quad = [1, 1, 1]
 72:   PetscInt, dimension(4), target      :: dofS3D = [0, 0, 0, 6]
 73:   PetscInt, dimension(4), target      :: dofUP1Tet = [3, 0, 0, 0]
 74:   PetscInt, dimension(4), target      :: dofAP1Tet = [1, 0, 0, 0]
 75:   PetscInt, dimension(4), target      :: dofUP2Tet = [3, 3, 0, 0]
 76:   PetscInt, dimension(4), target      :: dofAP2Tet = [1, 1, 0, 0]
 77:   PetscInt, dimension(4), target      :: dofUP1Hex = [3, 0, 0, 0]
 78:   PetscInt, dimension(4), target      :: dofAP1Hex = [1, 0, 0, 0]
 79:   PetscInt, dimension(4), target      :: dofUP2Hex = [3, 3, 3, 3]
 80:   PetscInt, dimension(4), target      :: dofAP2Hex = [1, 1, 1, 1]
 81:   PetscInt, dimension(:), pointer     :: dofU, dofA, dofS

 83:   type(tPetscSF)                      :: migrationSF
 84:   PetscPartitioner                    :: part
 85:   PetscLayout                         :: layout

 87:   type(tVec)                          :: tmpVec
 88:   PetscReal                           :: norm
 89:   PetscReal                           :: time

 91:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
 92:   if (ierr /= 0) then
 93:     print *, 'Unable to initialize PETSc'
 94:     stop
 95:   end if

 97:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
 98:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, numProc, ierr))
 99:   PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-i', ifilename, flg, ierr))
100:   PetscCheckA(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, 'missing input file name -i ')
101:   PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-o', ofilename, flg, ierr))
102:   PetscCheckA(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, 'missing output file name -o 
')
103:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-order', order, flg, ierr))
104:   if ((order > 2) .or. (order < 1)) then
105:     write (IOBuffer, '("Unsupported polynomial order ", I2, " not in [1,2]")') order
106:     SETERRA(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, IOBuffer)
107:   end if

109:   ! Read the mesh in any supported format
110:   PetscCallA(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, PETSC_NULL_CHARACTER, PETSC_TRUE, dm, ierr))
111:   PetscCallA(DMPlexDistributeSetDefault(dm, PETSC_FALSE, ierr))
112:   PetscCallA(DMSetFromOptions(dm, ierr))
113:   PetscCallA(DMGetDimension(dm, sdim, ierr))
114:   PetscCallA(DMViewFromOptions(dm, PETSC_NULL_OBJECT, '-dm_view', ierr))

116:   ! Create the exodus result file
117:   ! enable exodus debugging information
118:   PetscCallA(exopts(EXVRBS + EXDEBG, ierr))
119:   ! Create the exodus file
120:   PetscCallA(PetscViewerExodusIIOpen(PETSC_COMM_WORLD, ofilename, FILE_MODE_WRITE, viewer, ierr))
121:   ! The long way would be
122:   !
123:   ! PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD,viewer,ierr))
124:   ! PetscCallA(PetscViewerSetType(viewer,PETSCVIEWEREXODUSII,ierr))
125:   ! PetscCallA(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE,ierr))
126:   ! PetscCallA(PetscViewerFileSetName(viewer,ofilename,ierr))

128:   ! set the mesh order
129:   PetscCallA(PetscViewerExodusIISetOrder(viewer, order, ierr))
130:   PetscCallA(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD, ierr))
131:   !
132:   !    Notice how the exodus file is actually NOT open at this point (exoid is -1)
133:   !    Since we are overwriting the file (mode is FILE_MODE_WRITE), we are going to have to
134:   !    write the geometry (the DM), which can only be done on a brand new file.
135:   !

137:   ! Save the geometry to the file, erasing all previous content
138:   PetscCallA(DMView(dm, viewer, ierr))
139:   PetscCallA(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD, ierr))
140:   !
141:   !    Note how the exodus file is now open
142:   !
143:   ! 'Format' the exodus result file, i.e. allocate space for nodal and zonal variables
144:   select case (sdim)
145:   case (2)
146:     numNodalVar = 4
147:     numZonalVar = 3
148:     PetscCallA(PetscViewerExodusIISetZonalVariable(viewer, numZonalVar, ierr))
149:     PetscCallA(PetscViewerExodusIISetNodalVariable(viewer, numNodalVar, ierr))
150:     do i = 1, numZonalVar
151:       PetscCallA(PetscViewerExodusIISetZonalVariableName(viewer, i - 1, zonalVarName2D(i), ierr))
152:     end do
153:     do i = 1, numNodalVar
154:       PetscCallA(PetscViewerExodusIISetNodalVariableName(viewer, i - 1, nodalVarName2D(i), ierr))
155:     end do
156:   case (3)
157:     numNodalVar = 5
158:     numZonalVar = 6
159:     PetscCallA(PetscViewerExodusIISetZonalVariable(viewer, numZonalVar, ierr))
160:     PetscCallA(PetscViewerExodusIISetNodalVariable(viewer, numNodalVar, ierr))
161:     do i = 1, numZonalVar
162:       PetscCallA(PetscViewerExodusIISetZonalVariableName(viewer, i - 1, zonalVarName3D(i), ierr))
163:     end do
164:     do i = 1, numNodalVar
165:       PetscCallA(PetscViewerExodusIISetNodalVariableName(viewer, i - 1, nodalVarName3D(i), ierr))
166:     end do
167:   case default
168:     write (IOBuffer, '("No layout for dimension ",I2)') sdim
169:   end select
170:   PetscCallA(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD, ierr))

172:   !    An ExodusII truth table specifies which fields are saved at which time step
173:   !    It speeds up I/O but reserving space for fields in the file ahead of time.
174:   PetscCallA(PetscViewerExodusIIGetId(viewer, exoid, ierr))
175:   PetscCallA(exinq(exoid, EX_INQ_ELEM_BLK, numCS, PETSC_NULL_REAL, sjunk, ierr))
176:   allocate (truthtable(numCS, numZonalVar))
177:   truthtable = .true.
178:   PetscCallA(expvtt(exoid, numCS, numZonalVar, truthtable, ierr))
179:   deallocate (truthtable)

181:   !   Writing time step information in the file. Note that this is currently broken in the exodus library for netcdf4 (HDF5-based) file
182:   do step = 1, numstep
183:     time = real(step, kind=PETSC_REAL_KIND)
184:     PetscCallA(exptim(exoid, step, time, ierr))
185:   end do

187:   PetscCallA(PetscObjectGetComm(dm, comm, ierr))
188:   PetscCallA(PetscSectionCreate(comm, section, ierr))
189:   PetscCallA(PetscSectionSetNumFields(section, 3_PETSC_INT_KIND, ierr))
190:   PetscCallA(PetscSectionSetFieldName(section, fieldU, 'U', ierr))
191:   PetscCallA(PetscSectionSetFieldName(section, fieldA, 'Alpha', ierr))
192:   PetscCallA(PetscSectionSetFieldName(section, fieldS, 'Sigma', ierr))
193:   PetscCallA(DMPlexGetChart(dm, pStart, pEnd, ierr))
194:   PetscCallA(PetscSectionSetChart(section, pStart, pEnd, ierr))

196:   allocate (pStartDepth(sdim + 1))
197:   allocate (pEndDepth(sdim + 1))
198:   do d = 1, sdim + 1
199:     PetscCallA(DMPlexGetDepthStratum(dm, d - 1, pStartDepth(d), pEndDepth(d), ierr))
200:   end do

202:   ! Vector field U, Scalar field Alpha, Tensor field Sigma
203:   PetscCallA(PetscSectionSetFieldComponents(section, fieldU, sdim, ierr))
204:   PetscCallA(PetscSectionSetFieldComponents(section, fieldA, 1_PETSC_INT_KIND, ierr))
205:   PetscCallA(PetscSectionSetFieldComponents(section, fieldS, sdim*(sdim + 1)/2, ierr))

207:   ! Going through cell sets then cells, and setting up storage for the sections
208:   PetscCallA(DMGetLabelSize(dm, 'Cell Sets', numCS, ierr))
209:   PetscCallA(DMGetLabelIdIS(dm, 'Cell Sets', csIS, ierr))
210:   PetscCallA(ISGetIndices(csIS, csID, ierr))
211:   do set = 1, numCS
212:         !!! This is pointless but will avoid a warning with gfortran and -Werror=maybe-uninitialized
213:     nullify (dofA)
214:     nullify (dofU)
215:     nullify (dofS)
216:     PetscCallA(DMGetStratumSize(dm, 'Cell Sets', csID(set), numCells, ierr))
217:     PetscCallA(DMGetStratumIS(dm, 'Cell Sets', csID(set), cellIS, ierr))
218:     if (numCells > 0) then
219:       select case (sdim)
220:       case (2)
221:         dofs => dofS2D
222:       case (3)
223:         dofs => dofS3D
224:       case default
225:         write (IOBuffer, '("No layout for dimension ",I2)') sdim
226:         SETERRA(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, IOBuffer)
227:       end select ! sdim

229:       ! Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes
230:       ! It will not be enough to identify more exotic elements like pyramid or prisms...  */
231:       PetscCallA(ISGetIndices(cellIS, cellID, ierr))
232:       nullify (closureA)
233:       PetscCallA(DMPlexGetTransitiveClosure(dm, cellID(1), PETSC_TRUE, PETSC_NULL_INTEGER, closureA, ierr))
234:       select case (size(closureA)/2)
235:       case (7) ! Tri
236:         if (order == 1) then
237:           dofU => dofUP1Tri
238:           dofA => dofAP1Tri
239:         else
240:           dofU => dofUP2Tri
241:           dofA => dofAP2Tri
242:         end if
243:       case (9) ! Quad
244:         if (order == 1) then
245:           dofU => dofUP1Quad
246:           dofA => dofAP1Quad
247:         else
248:           dofU => dofUP2Quad
249:           dofA => dofAP2Quad
250:         end if
251:       case (15) ! Tet
252:         if (order == 1) then
253:           dofU => dofUP1Tet
254:           dofA => dofAP1Tet
255:         else
256:           dofU => dofUP2Tet
257:           dofA => dofAP2Tet
258:         end if
259:       case (27) ! Hex
260:         if (order == 1) then
261:           dofU => dofUP1Hex
262:           dofA => dofAP1Hex
263:         else
264:           dofU => dofUP2Hex
265:           dofA => dofAP2Hex
266:         end if
267:       case default
268:         write (IOBuffer, '("Unknown element with closure size ",I2)') size(closureA)/2
269:         SETERRA(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, IOBuffer)
270:       end select
271:       PetscCallA(DMPlexRestoreTransitiveClosure(dm, cellID(1), PETSC_TRUE, PETSC_NULL_INTEGER, closureA, ierr))
272:       do cell = 1, numCells!
273:         nullify (closure)
274:         PetscCallA(DMPlexGetTransitiveClosure(dm, cellID(cell), PETSC_TRUE, PETSC_NULL_INTEGER, closure, ierr))
275:         do p = 1, size(closure), 2
276:           ! find the depth of p
277:           do d = 1, sdim + 1
278:             if ((closure(p) >= pStartDepth(d)) .and. (closure(p) < pEndDepth(d))) then
279:               PetscCallA(PetscSectionSetDof(section, closure(p), dofU(d) + dofA(d) + dofS(d), ierr))
280:               PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldU, dofU(d), ierr))
281:               PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldA, dofA(d), ierr))
282:               PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldS, dofS(d), ierr))
283:             end if ! closure(p)
284:           end do ! d
285:         end do ! p
286:         PetscCallA(DMPlexRestoreTransitiveClosure(dm, cellID(cell), PETSC_TRUE, PETSC_NULL_INTEGER, closure, ierr))
287:       end do ! cell
288:       PetscCallA(ISRestoreIndices(cellIS, cellID, ierr))
289:       PetscCallA(ISDestroy(cellIS, ierr))
290:     end if ! numCells
291:   end do ! set
292:   PetscCallA(ISRestoreIndices(csIS, csID, ierr))
293:   PetscCallA(ISDestroy(csIS, ierr))
294:   PetscCallA(PetscSectionSetUp(section, ierr))
295:   PetscCallA(DMSetLocalSection(dm, section, ierr))
296:   PetscCallA(PetscObjectViewFromOptions(section, PETSC_NULL_OBJECT, '-dm_section_view', ierr))
297:   PetscCallA(PetscSectionGetValueLayout(PETSC_COMM_WORLD, section, layout, ierr))
298:   PetscCallA(PetscLayoutDestroy(layout, ierr))
299:   PetscCallA(PetscSectionDestroy(section, ierr))

301:   PetscCallA(DMSetUseNatural(dm, PETSC_TRUE, ierr))
302:   PetscCallA(DMPlexGetPartitioner(dm, part, ierr))
303:   PetscCallA(PetscPartitionerSetFromOptions(part, ierr))
304:   PetscCallA(DMPlexDistribute(dm, 0_PETSC_INT_KIND, migrationSF, pdm, ierr))

306:   if (numProc > 1) then
307:     PetscCallA(DMPlexSetMigrationSF(pdm, migrationSF, ierr))
308:     PetscCallA(PetscSFDestroy(migrationSF, ierr))
309:   else
310:     pdm = dm
311:   end if
312:   PetscCallA(DMViewFromOptions(pdm, PETSC_NULL_OBJECT, '-dm_view', ierr))

314:   ! Get DM and IS for each field of dm
315:   PetscCallA(DMCreateSubDM(pdm, 1_PETSC_INT_KIND, [fieldU], isU, dmU, ierr))
316:   PetscCallA(DMCreateSubDM(pdm, 1_PETSC_INT_KIND, [fieldA], isA, dmA, ierr))
317:   PetscCallA(DMCreateSubDM(pdm, 1_PETSC_INT_KIND, [fieldS], isS, dmS, ierr))
318:   PetscCallA(DMCreateSubDM(pdm, 2_PETSC_INT_KIND, fieldUA, isUA, dmUA, ierr))

320:   !Create the exodus result file
321:   allocate (dmList(2))
322:   dmList(1) = dmU
323:   dmList(2) = dmA
324:   PetscCallA(DMCreateSuperDM(dmList, 2_PETSC_INT_KIND, PETSC_NULL_IS_POINTER, dmUA2, ierr))
325:   deallocate (dmList)

327:   PetscCallA(DMGetGlobalVector(pdm, X, ierr))
328:   PetscCallA(DMGetGlobalVector(dmU, U, ierr))
329:   PetscCallA(DMGetGlobalVector(dmA, A, ierr))
330:   PetscCallA(DMGetGlobalVector(dmS, S, ierr))
331:   PetscCallA(DMGetGlobalVector(dmUA, UA, ierr))
332:   PetscCallA(DMGetGlobalVector(dmUA2, UA2, ierr))

334:   PetscCallA(PetscObjectSetName(U, 'U', ierr))
335:   PetscCallA(PetscObjectSetName(A, 'Alpha', ierr))
336:   PetscCallA(PetscObjectSetName(S, 'Sigma', ierr))
337:   PetscCallA(PetscObjectSetName(UA, 'UAlpha', ierr))
338:   PetscCallA(PetscObjectSetName(UA2, 'UAlpha2', ierr))
339:   PetscCallA(VecSet(X, -111.0_PETSC_REAL_KIND, ierr))

341:   ! 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 */
342:   PetscCallA(DMGetLocalSection(dmUA, sectionUA, ierr))
343:   PetscCallA(DMGetLocalVector(dmUA, UALoc, ierr))
344:   PetscCallA(VecGetArray(UALoc, cval, ierr))
345:   PetscCallA(DMGetCoordinateSection(dmUA, coordSection, ierr))
346:   PetscCallA(DMGetCoordinatesLocal(dmUA, coord, ierr))
347:   PetscCallA(DMPlexGetChart(dmUA, pStart, pEnd, ierr))

349:   do p = pStart, pEnd - 1
350:     PetscCallA(PetscSectionGetDof(sectionUA, p, dofUA, ierr))
351:     if (dofUA > 0) then
352:       PetscCallA(PetscSectionGetOffset(sectionUA, p, offUA, ierr))
353:       PetscCallA(DMPlexVecGetClosure(dmUA, coordSection, coord, p, PETSC_NULL_INTEGER, xyz, ierr))
354:       closureSize = size(xyz)
355:       do k = 1, sdim
356:         do j = 0, closureSize - 1, sdim
357:           cval(offUA + k) = cval(offUA + k) + xyz(j/sdim + k)
358:         end do
359:         cval(offUA + k) = cval(offUA + k)*sdim/closureSize
360:         cval(offUA + sdim + 1) = cval(offUA + sdim + 1) + cval(offUA + k)**2
361:       end do
362:       PetscCallA(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, PETSC_NULL_INTEGER, xyz, ierr))
363:     end if
364:   end do

366:   PetscCallA(VecRestoreArray(UALoc, cval, ierr))
367:   PetscCallA(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA, ierr))
368:   PetscCallA(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA, ierr))
369:   PetscCallA(DMRestoreLocalVector(dmUA, UALoc, ierr))

371:   !Update X
372:   PetscCallA(VecISCopy(X, isUA, SCATTER_FORWARD, UA, ierr))
373:   ! Restrict to U and Alpha
374:   PetscCallA(VecISCopy(X, isU, SCATTER_REVERSE, U, ierr))
375:   PetscCallA(VecISCopy(X, isA, SCATTER_REVERSE, A, ierr))
376:   PetscCallA(VecViewFromOptions(UA, PETSC_NULL_OBJECT, '-ua_vec_view', ierr))
377:   PetscCallA(VecViewFromOptions(U, PETSC_NULL_OBJECT, '-u_vec_view', ierr))
378:   PetscCallA(VecViewFromOptions(A, PETSC_NULL_OBJECT, '-a_vec_view', ierr))
379:   ! restrict to UA2
380:   PetscCallA(VecISCopy(X, isUA, SCATTER_REVERSE, UA2, ierr))
381:   PetscCallA(VecViewFromOptions(UA2, PETSC_NULL_OBJECT, '-ua2_vec_view', ierr))

383:   ! Getting Natural Vec
384:   PetscCallA(DMSetOutputSequenceNumber(dmU, 0_PETSC_INT_KIND, time, ierr))
385:   PetscCallA(DMSetOutputSequenceNumber(dmA, 0_PETSC_INT_KIND, time, ierr))

387:   PetscCallA(VecView(U, viewer, ierr))
388:   PetscCallA(VecView(A, viewer, ierr))

390:   !Saving U and Alpha in one shot.
391:   !For this, we need to cheat and change the Vec's name
392:   !Note that in the end we write variables one component at a time,
393:   !so that there is no real value in doing this
394:   PetscCallA(DMSetOutputSequenceNumber(dmUA, 1_PETSC_INT_KIND, time, ierr))
395:   PetscCallA(DMGetGlobalVector(dmUA, tmpVec, ierr))
396:   PetscCallA(VecCopy(UA, tmpVec, ierr))
397:   PetscCallA(PetscObjectSetName(tmpVec, 'U', ierr))
398:   PetscCallA(VecView(tmpVec, viewer, ierr))

400:   ! Reading nodal variables in Exodus file
401:   PetscCallA(VecSet(tmpVec, -1000.0_PETSC_REAL_KIND, ierr))
402:   PetscCallA(VecLoad(tmpVec, viewer, ierr))
403:   PetscCallA(VecAXPY(UA, -1.0_PETSC_REAL_KIND, tmpVec, ierr))
404:   PetscCallA(VecNorm(UA, NORM_INFINITY, norm, ierr))
405:   if (norm > PETSC_SQRT_MACHINE_EPSILON) then
406:     write (IOBuffer, '("UAlpha ||Vin - Vout|| = ",ES12.5)') norm
407:   end if
408:   PetscCallA(DMRestoreGlobalVector(dmUA, tmpVec, ierr))

410:   ! same thing with the UA2 Vec obtained from the superDM
411:   PetscCallA(DMGetGlobalVector(dmUA2, tmpVec, ierr))
412:   PetscCallA(VecCopy(UA2, tmpVec, ierr))
413:   PetscCallA(PetscObjectSetName(tmpVec, 'U', ierr))
414:   PetscCallA(DMSetOutputSequenceNumber(dmUA2, 2_PETSC_INT_KIND, time, ierr))
415:   PetscCallA(VecView(tmpVec, viewer, ierr))

417:   ! Reading nodal variables in Exodus file
418:   PetscCallA(VecSet(tmpVec, -1000.0_PETSC_REAL_KIND, ierr))
419:   PetscCallA(VecLoad(tmpVec, viewer, ierr))
420:   PetscCallA(VecAXPY(UA2, -1.0_PETSC_REAL_KIND, tmpVec, ierr))
421:   PetscCallA(VecNorm(UA2, NORM_INFINITY, norm, ierr))
422:   if (norm > PETSC_SQRT_MACHINE_EPSILON) then
423:     write (IOBuffer, '("UAlpha2 ||Vin - Vout|| = ",ES12.5)') norm
424:   end if
425:   PetscCallA(DMRestoreGlobalVector(dmUA2, tmpVec, ierr))

427:   ! Building and saving Sigma
428:   !   We set sigma_0 = rank (to see partitioning)
429:   !          sigma_1 = cell set ID
430:   !          sigma_2 = x_coordinate of the cell center of mass
431:   PetscCallA(DMGetCoordinateSection(dmS, coordSection, ierr))
432:   PetscCallA(DMGetCoordinatesLocal(dmS, coord, ierr))
433:   PetscCallA(DMGetLabelIdIS(dmS, 'Cell Sets', csIS, ierr))
434:   PetscCallA(DMGetLabelSize(dmS, 'Cell Sets', numCS, ierr))
435:   PetscCallA(ISGetIndices(csIS, csID, ierr))

437:   do set = 1, numCS
438:     PetscCallA(DMGetStratumIS(dmS, 'Cell Sets', csID(set), cellIS, ierr))
439:     PetscCallA(ISGetIndices(cellIS, cellID, ierr))
440:     PetscCallA(ISGetSize(cellIS, numCells, ierr))
441:     do cell = 1, numCells
442:       PetscCallA(DMPlexVecGetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), PETSC_NULL_INTEGER, cval, ierr))
443:       PetscCallA(DMPlexVecGetClosure(dmS, coordSection, coord, cellID(cell), PETSC_NULL_INTEGER, xyz, ierr))
444:       cval(1) = rank
445:       cval(2) = csID(set)
446:       cval(3) = 0.0_PETSC_REAL_KIND
447:       do c = 1, size(xyz), sdim
448:         cval(3) = cval(3) + xyz(c)
449:       end do
450:       cval(3) = cval(3)*sdim/size(xyz)
451:       PetscCallA(DMPlexVecSetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval, INSERT_ALL_VALUES, ierr))
452:       PetscCallA(DMPlexVecRestoreClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), PETSC_NULL_INTEGER, cval, ierr))
453:       PetscCallA(DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID(cell), PETSC_NULL_INTEGER, xyz, ierr))
454:     end do
455:     PetscCallA(ISRestoreIndices(cellIS, cellID, ierr))
456:     PetscCallA(ISDestroy(cellIS, ierr))
457:   end do
458:   PetscCallA(ISRestoreIndices(csIS, csID, ierr))
459:   PetscCallA(ISDestroy(csIS, ierr))
460:   PetscCallA(VecViewFromOptions(S, PETSC_NULL_OBJECT, '-s_vec_view', ierr))

462:   ! Writing zonal variables in Exodus file
463:   PetscCallA(DMSetOutputSequenceNumber(dmS, 0_PETSC_INT_KIND, time, ierr))
464:   PetscCallA(VecView(S, viewer, ierr))

466:   ! Reading zonal variables in Exodus file */
467:   PetscCallA(DMGetGlobalVector(dmS, tmpVec, ierr))
468:   PetscCallA(VecSet(tmpVec, -1000.0_PETSC_REAL_KIND, ierr))
469:   PetscCallA(PetscObjectSetName(tmpVec, 'Sigma', ierr))
470:   PetscCallA(VecLoad(tmpVec, viewer, ierr))
471:   PetscCallA(VecAXPY(S, -1.0_PETSC_REAL_KIND, tmpVec, ierr))
472:   PetscCallA(VecNorm(S, NORM_INFINITY, norm, ierr))
473:   if (norm > PETSC_SQRT_MACHINE_EPSILON) then
474:     write (IOBuffer, '("Sigma ||Vin - Vout|| = ",ES12.5)') norm
475:   end if
476:   PetscCallA(DMRestoreGlobalVector(dmS, tmpVec, ierr))

478:   PetscCallA(DMRestoreGlobalVector(dmUA2, UA2, ierr))
479:   PetscCallA(DMRestoreGlobalVector(dmUA, UA, ierr))
480:   PetscCallA(DMRestoreGlobalVector(dmS, S, ierr))
481:   PetscCallA(DMRestoreGlobalVector(dmA, A, ierr))
482:   PetscCallA(DMRestoreGlobalVector(dmU, U, ierr))
483:   PetscCallA(DMRestoreGlobalVector(pdm, X, ierr))
484:   PetscCallA(DMDestroy(dmU, ierr))
485:   PetscCallA(ISDestroy(isU, ierr))
486:   PetscCallA(DMDestroy(dmA, ierr))
487:   PetscCallA(ISDestroy(isA, ierr))
488:   PetscCallA(DMDestroy(dmS, ierr))
489:   PetscCallA(ISDestroy(isS, ierr))
490:   PetscCallA(DMDestroy(dmUA, ierr))
491:   PetscCallA(ISDestroy(isUA, ierr))
492:   PetscCallA(DMDestroy(dmUA2, ierr))
493:   PetscCallA(DMDestroy(pdm, ierr))
494:   if (numProc > 1) then
495:     PetscCallA(DMDestroy(dm, ierr))
496:   end if

498:   deallocate (pStartDepth)
499:   deallocate (pEndDepth)

501:   PetscCallA(PetscViewerDestroy(viewer, ierr))
502:   PetscCallA(PetscFinalize(ierr))
503: end program ex26f90

505: ! /*TEST
506: !
507: ! build:
508: !   requires: exodusii pnetcdf !complex
509: !   # 2D seq
510: ! test:
511: !   suffix: 0
512: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
513: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
514: ! test:
515: !   suffix: 1
516: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
517: !
518: ! test:
519: !   suffix: 2
520: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
521: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
522: ! test:
523: !   suffix: 3
524: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
525: ! test:
526: !   suffix: 4
527: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
528: ! test:
529: !   suffix: 5
530: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
531: !   # 2D par
532: ! test:
533: !   suffix: 6
534: !   nsize: 2
535: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
536: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
537: ! test:
538: !   suffix: 7
539: !   nsize: 2
540: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
541: ! test:
542: !   suffix: 8
543: !   nsize: 2
544: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
545: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name
546: ! test:
547: !   suffix: 9
548: !   nsize: 2
549: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
550: ! test:
551: !   suffix: 10
552: !   nsize: 2
553: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
554: ! test:
555: !   # Something is now broken with parallel read/write for whatever shape H is
556: !   TODO: broken
557: !   suffix: 11
558: !   nsize: 2
559: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2

561: !   #3d seq
562: ! test:
563: !   suffix: 12
564: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
565: ! test:
566: !   suffix: 13
567: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
568: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
569: ! test:
570: !   suffix: 14
571: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
572: ! test:
573: !   suffix: 15
574: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
575: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
576: !   #3d par
577: ! test:
578: !   suffix: 16
579: !   nsize: 2
580: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
581: ! test:
582: !   suffix: 17
583: !   nsize: 2
584: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 1
585: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints
586: ! test:
587: !   suffix: 18
588: !   nsize: 2
589: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
590: ! test:
591: !   suffix: 19
592: !   nsize: 2
593: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -dm_section_view -petscpartitioner_type simple -order 2
594: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints

596: ! TEST*/