1: #include "petsc/finclude/petsc.h"
  2: program ex62f90
  3:   use petsc
  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, rootSection, leafSection
 12:   PetscInt, parameter                :: fieldU = 0, fieldS = 1, fieldA = 2
 13:   PetscInt, 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:   PetscInt                           :: sdim, d, pStart, pEnd, p, numCS, set, i, j
 21:   PetscMPIInt                        :: rank, numProc
 22:   PetscBool                          :: flg
 23:   PetscErrorCode                     :: ierr
 24:   MPIU_Comm                          :: comm
 25:   type(tPetscViewer)                 :: viewer

 27:   character(len=MXSTLN)              :: sJunk
 28:   PetscInt, parameter                :: numstep = 3
 29:   PetscInt                           :: step
 30:   PetscInt                           :: numNodalVar, numZonalVar
 31:   character(len=MXSTLN)              :: nodalVarName(4)
 32:   character(len=MXSTLN)              :: zonalVarName(6)
 33:   logical, dimension(:, :), pointer  :: truthtable

 35:   type(tIS)                          :: cellIS
 36:   PetscInt, dimension(:), pointer    :: cellID
 37:   PetscInt                           :: numCells, cell, closureSize
 38:   PetscInt, dimension(:), pointer    :: closureA, closure

 40:   type(tPetscSection)                :: sectionUA, coordSection
 41:   type(tVec)                         :: UALoc, coord
 42:   PetscScalar, dimension(:), pointer :: cval, xyz
 43:   PetscInt                           :: dofUA, offUA, c

 45:   ! dof layout ordered by increasing height in the DAG: cell, face, edge, vertex
 46:   PetscInt, dimension(3), target     :: dofS2D = [0, 0, 3]
 47:   PetscInt, dimension(3), target     :: dofUP1Tri = [2, 0, 0]
 48:   PetscInt, dimension(3), target     :: dofAP1Tri = [1, 0, 0]
 49:   PetscInt, dimension(3), target     :: dofUP2Tri = [2, 2, 0]
 50:   PetscInt, dimension(3), target     :: dofAP2Tri = [1, 1, 0]
 51:   PetscInt, dimension(3), target     :: dofUP1Quad = [2, 0, 0]
 52:   PetscInt, dimension(3), target     :: dofAP1Quad = [1, 0, 0]
 53:   PetscInt, dimension(3), target     :: dofUP2Quad = [2, 2, 2]
 54:   PetscInt, dimension(3), target     :: dofAP2Quad = [1, 1, 1]
 55:   PetscInt, dimension(4), target     :: dofS3D = [0, 0, 0, 6]
 56:   PetscInt, dimension(4), target     :: dofUP1Tet = [3, 0, 0, 0]
 57:   PetscInt, dimension(4), target     :: dofAP1Tet = [1, 0, 0, 0]
 58:   PetscInt, dimension(4), target     :: dofUP2Tet = [3, 3, 0, 0]
 59:   PetscInt, dimension(4), target     :: dofAP2Tet = [1, 1, 0, 0]
 60:   PetscInt, dimension(4), target     :: dofUP1Hex = [3, 0, 0, 0]
 61:   PetscInt, dimension(4), target     :: dofAP1Hex = [1, 0, 0, 0]
 62:   PetscInt, dimension(4), target     :: dofUP2Hex = [3, 3, 3, 3]
 63:   PetscInt, dimension(4), target     :: dofAP2Hex = [1, 1, 1, 1]
 64:   PetscInt, dimension(:), pointer    :: dofU, dofA, dofS

 66:   type(tPetscSF)                      :: migrationSF, natSF, natPointSF, natPointSFInv
 67:   PetscPartitioner                    :: part

 69:   type(tVec)                          :: tmpVec
 70:   PetscReal                           :: norm
 71:   PetscReal                           :: time
 72:   PetscInt, pointer                   :: remoteOffsets(:)

 74:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
 75:   if (ierr /= 0) then
 76:     print *, 'Unable to initialize PETSc'
 77:     stop
 78:   end if

 80:   PetscCallA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
 81:   PetscCallA(MPI_Comm_size(PETSC_COMM_WORLD, numProc, ierr))
 82:   PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-i', ifilename, flg, ierr))
 83:   PetscCheckA(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, 'missing input file name -i ')
 84:   PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-o', ofilename, flg, ierr))
 85:   PetscCheckA(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, 'missing output file name -o 
')
 86:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-order', order, flg, ierr))
 87:   if ((order > 2) .or. (order < 1)) then
 88:     write (IOBuffer, '("Unsupported polynomial order ", I2, " not in [1,2]")') order
 89:     stop
 90:   end if

 92:   ! Read the mesh in any supported format
 93:   PetscCallA(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, PETSC_NULL_CHARACTER, PETSC_TRUE, dm, ierr))
 94:   PetscCallA(DMPlexDistributeSetDefault(dm, PETSC_FALSE, ierr))
 95:   PetscCallA(DMSetFromOptions(dm, ierr))
 96:   PetscCallA(DMGetDimension(dm, sdim, ierr))
 97:   PetscCallA(DMViewFromOptions(dm, PETSC_NULL_OBJECT, '-dm_view', ierr))

 99:   ! Create the exodus result file

101:   ! enable exodus debugging information
102:   PetscCallA(exopts(EXVRBS + EXDEBG, ierr))
103:   ! Create the exodus file
104:   PetscCallA(PetscViewerExodusIIOpen(PETSC_COMM_WORLD, ofilename, FILE_MODE_WRITE, viewer, ierr))
105:   ! The long way would be
106:   !
107:   ! PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD,viewer,ierr))
108:   ! PetscCallA(PetscViewerSetType(viewer,PETSCVIEWEREXODUSII,ierr))
109:   ! PetscCallA(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE,ierr))
110:   ! PetscCallA(PetscViewerFileSetName(viewer,ofilename,ierr))

112:   ! set the mesh order
113:   PetscCallA(PetscViewerExodusIISetOrder(viewer, order, ierr))
114:   PetscCallA(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD, ierr))
115:   !
116:   !    Notice how the exodus file is actually NOT open at this point (exoid is -1)
117:   !    Since we are overwriting the file (mode is FILE_MODE_WRITE), we are going to have to
118:   !    write the geometry (the DM), which can only be done on a brand new file.
119:   !

121:   ! Save the geometry to the file, erasing all previous content
122:   PetscCallA(DMView(dm, viewer, ierr))
123:   PetscCallA(PetscViewerView(viewer, PETSC_VIEWER_STDOUT_WORLD, ierr))
124:   !
125:   !    Note how the exodus file is now open
126:   !
127:   ! 'Format' the exodus result file, i.e. allocate space for nodal and zonal variables
128:   select case (sdim)
129:   case (2)
130:     numNodalVar = 3
131:     nodalVarName(1:numNodalVar) = ['U_x  ', 'U_y  ', 'Alpha']
132:     numZonalVar = 3
133:     zonalVarName(1:numZonalVar) = ['Sigma_11', 'Sigma_22', 'Sigma_12']
134:   case (3)
135:     numNodalVar = 4
136:     nodalVarName(1:numNodalVar) = ['U_x  ', 'U_y  ', 'U_z  ', 'Alpha']
137:     numZonalVar = 6
138:     zonalVarName(1:numZonalVar) = ['Sigma_11', 'Sigma_22', 'Sigma_33', 'Sigma_23', 'Sigma_13', 'Sigma_12']
139:   case default
140:     write (IOBuffer, '("No layout for dimension ",I2)') sdim
141:   end select
142:   PetscCallA(PetscViewerExodusIIGetId(viewer, exoid, ierr))
143:   PetscCallA(expvp(exoid, 'E', numZonalVar, ierr))
144:   PetscCallA(expvan(exoid, 'E', numZonalVar, zonalVarName, ierr))
145:   PetscCallA(expvp(exoid, 'N', numNodalVar, ierr))
146:   PetscCallA(expvan(exoid, 'N', numNodalVar, nodalVarName, ierr))
147:   PetscCallA(exinq(exoid, EX_INQ_ELEM_BLK, numCS, PETSC_NULL_REAL, sjunk, ierr))

149:   !    An exodusII truth table specifies which fields are saved at which time step
150:   !    It speeds up I/O but reserving space for fields in the file ahead of time.
151:   allocate (truthtable(numCS, numZonalVar))
152:   truthtable = .true.
153:   PetscCallA(expvtt(exoid, numCS, numZonalVar, truthtable, ierr))
154:   deallocate (truthtable)

156:   !   Writing time step information in the file. Note that this is currently broken in the exodus library for netcdf4 (HDF5-based) files
157:   do step = 1, numstep
158:     time = real(step, kind=PETSC_REAL_KIND)
159:     PetscCallA(exptim(exoid, step, time, ierr))
160:   end do

162:   PetscCallA(DMSetUseNatural(dm, PETSC_TRUE, ierr))
163:   PetscCallA(DMPlexGetPartitioner(dm, part, ierr))
164:   PetscCallA(PetscPartitionerSetFromOptions(part, ierr))
165:   PetscCallA(DMPlexDistribute(dm, 0_PETSC_INT_KIND, migrationSF, pdm, ierr))

167:   if (numProc > 1) then
168:     PetscCallA(DMPlexSetMigrationSF(pdm, migrationSF, ierr))
169:     PetscCallA(PetscSFDestroy(migrationSF, ierr))
170:   else
171:     pdm = dm
172:   end if
173:   PetscCallA(DMViewFromOptions(pdm, PETSC_NULL_OBJECT, '-dm_view', ierr))

175:   PetscCallA(PetscObjectGetComm(pdm, comm, ierr))
176:   PetscCallA(PetscSectionCreate(comm, section, ierr))
177:   PetscCallA(PetscSectionSetNumFields(section, 3_PETSC_INT_KIND, ierr))
178:   PetscCallA(PetscSectionSetFieldName(section, fieldU, 'U', ierr))
179:   PetscCallA(PetscSectionSetFieldName(section, fieldA, 'Alpha', ierr))
180:   PetscCallA(PetscSectionSetFieldName(section, fieldS, 'Sigma', ierr))
181:   PetscCallA(DMPlexGetChart(pdm, pStart, pEnd, ierr))
182:   PetscCallA(PetscSectionSetChart(section, pStart, pEnd, ierr))

184:   allocate (pStartDepth(sdim + 1))
185:   allocate (pEndDepth(sdim + 1))
186:   do d = 1, sdim + 1
187:     PetscCallA(DMPlexGetDepthStratum(pdm, d - 1, pStartDepth(d), pEndDepth(d), ierr))
188:   end do

190:   ! Vector field U, Scalar field Alpha, Tensor field Sigma
191:   PetscCallA(PetscSectionSetFieldComponents(section, fieldU, sdim, ierr))
192:   PetscCallA(PetscSectionSetFieldComponents(section, fieldA, 1_PETSC_INT_KIND, ierr))
193:   PetscCallA(PetscSectionSetFieldComponents(section, fieldS, sdim*(sdim + 1)/2, ierr))

195:   ! Going through cell sets then cells, and setting up storage for the sections
196:   PetscCallA(DMGetLabelSize(pdm, 'Cell Sets', numCS, ierr))
197:   PetscCallA(DMGetLabelIdIS(pdm, 'Cell Sets', csIS, ierr))
198:   PetscCallA(ISGetIndices(csIS, csID, ierr))
199:   do set = 1, numCS
200:     PetscCallA(DMGetStratumSize(pdm, 'Cell Sets', csID(set), numCells, ierr))
201:     PetscCallA(DMGetStratumIS(pdm, 'Cell Sets', csID(set), cellIS, ierr))
202:     if (numCells > 0) then
203:       select case (sdim)
204:       case (2)
205:         dofs => dofS2D
206:       case (3)
207:         dofs => dofS3D
208:       case default
209:         write (IOBuffer, '("No layout for dimension ",I2)') sdim
210:         stop
211:       end select ! sdim

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

283:   ! Creating section on the sequential DM + creating the GlobalToNatural SF
284:   if (numProc > 1) then
285:     PetscCallA(DMPlexGetMigrationSF(pdm, natPointSF, ierr))
286:     PetscCallA(DMGetLocalSection(pdm, rootSection, ierr))
287:     PetscCallA(PetscSFCreateInverseSF(natPointSF, natPointSFInv, ierr))
288:     PetscCallA(PetscSectionCreate(PETSC_COMM_WORLD, leafSection, ierr))
289:     PetscCallA(PetscSFDistributeSection(natPointSFInv, rootSection, remoteOffsets, leafSection, ierr))
290:     PetscCallA(PetscSFDestroyRemoteOffsets(remoteOffsets, ierr))
291:     PetscCallA(DMSetLocalSection(dm, leafSection, ierr))
292:     PetscCallA(DMPlexCreateGlobalToNaturalSF(pdm, leafSection, natPointSF, natSF, ierr))
293:     PetscCallA(PetscSFDestroy(natPointSFInv, ierr))
294:     PetscCallA(PetscSectionDestroy(leafSection, ierr))
295:     PetscCallA(DMSetNaturalSF(pdm, natSF, ierr))
296:     PetscCallA(PetscObjectDereference(natSF, ierr))
297:   end if

299:   ! Get DM and IS for each field of dm
300:   PetscCallA(DMCreateSubDM(pdm, 1_PETSC_INT_KIND, [fieldU], isU, dmU, ierr))
301:   PetscCallA(DMCreateSubDM(pdm, 1_PETSC_INT_KIND, [fieldA], isA, dmA, ierr))
302:   PetscCallA(DMCreateSubDM(pdm, 1_PETSC_INT_KIND, [fieldS], isS, dmS, ierr))
303:   PetscCallA(DMCreateSubDM(pdm, 2_PETSC_INT_KIND, fieldUA, isUA, dmUA, ierr))

305:   !Create the exodus result file
306:   allocate (dmList(2))
307:   dmList(1) = dmU
308:   dmList(2) = dmA
309:   PetscCallA(DMCreateSuperDM(dmList, 2_PETSC_INT_KIND, PETSC_NULL_IS_POINTER, dmUA2, ierr))
310:   deallocate (dmList)

312:   PetscCallA(DMGetGlobalVector(pdm, X, ierr))
313:   PetscCallA(DMGetGlobalVector(dmU, U, ierr))
314:   PetscCallA(DMGetGlobalVector(dmA, A, ierr))
315:   PetscCallA(DMGetGlobalVector(dmS, S, ierr))
316:   PetscCallA(DMGetGlobalVector(dmUA, UA, ierr))
317:   PetscCallA(DMGetGlobalVector(dmUA2, UA2, ierr))

319:   PetscCallA(PetscObjectSetName(U, 'U', ierr))
320:   PetscCallA(PetscObjectSetName(A, 'Alpha', ierr))
321:   PetscCallA(PetscObjectSetName(S, 'Sigma', ierr))
322:   PetscCallA(PetscObjectSetName(UA, 'UAlpha', ierr))
323:   PetscCallA(PetscObjectSetName(UA2, 'UAlpha2', ierr))
324:   PetscCallA(VecSet(X, -111.0_PETSC_REAL_KIND, ierr))

326:   ! 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 */
327:   PetscCallA(DMGetLocalSection(dmUA, sectionUA, ierr))
328:   PetscCallA(DMGetLocalVector(dmUA, UALoc, ierr))
329:   PetscCallA(VecGetArray(UALoc, cval, ierr))
330:   PetscCallA(DMGetCoordinateSection(dmUA, coordSection, ierr))
331:   PetscCallA(DMGetCoordinatesLocal(dmUA, coord, ierr))
332:   PetscCallA(DMPlexGetChart(dmUA, pStart, pEnd, ierr))

334:   do p = pStart, pEnd - 1
335:     PetscCallA(PetscSectionGetDof(sectionUA, p, dofUA, ierr))
336:     if (dofUA > 0) then
337:       PetscCallA(PetscSectionGetOffset(sectionUA, p, offUA, ierr))
338:       PetscCallA(DMPlexVecGetClosure(dmUA, coordSection, coord, p, PETSC_NULL_INTEGER, xyz, ierr))
339:       closureSize = size(xyz)
340:       do i = 1, sdim
341:         do j = 0, closureSize - 1, sdim
342:           cval(offUA + i) = cval(offUA + i) + xyz(j/sdim + i)
343:         end do
344:         cval(offUA + i) = cval(offUA + i)*sdim/closureSize
345:         cval(offUA + sdim + 1) = cval(offUA + sdim + 1) + cval(offUA + i)**2
346:       end do
347:       PetscCallA(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, PETSC_NULL_INTEGER, xyz, ierr))
348:     end if
349:   end do

351:   PetscCallA(VecRestoreArray(UALoc, cval, ierr))
352:   PetscCallA(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA, ierr))
353:   PetscCallA(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA, ierr))
354:   PetscCallA(DMRestoreLocalVector(dmUA, UALoc, ierr))

356:   !Update X
357:   PetscCallA(VecISCopy(X, isUA, SCATTER_FORWARD, UA, ierr))
358:   ! Restrict to U and Alpha
359:   PetscCallA(VecISCopy(X, isU, SCATTER_REVERSE, U, ierr))
360:   PetscCallA(VecISCopy(X, isA, SCATTER_REVERSE, A, ierr))
361:   PetscCallA(VecViewFromOptions(UA, PETSC_NULL_OBJECT, '-ua_vec_view', ierr))
362:   PetscCallA(VecViewFromOptions(U, PETSC_NULL_OBJECT, '-u_vec_view', ierr))
363:   PetscCallA(VecViewFromOptions(A, PETSC_NULL_OBJECT, '-a_vec_view', ierr))
364:   ! restrict to UA2
365:   PetscCallA(VecISCopy(X, isUA, SCATTER_REVERSE, UA2, ierr))
366:   PetscCallA(VecViewFromOptions(UA2, PETSC_NULL_OBJECT, '-ua2_vec_view', ierr))

368:   ! Getting Natural Vec
369:   PetscCallA(DMSetOutputSequenceNumber(dmU, 0_PETSC_INT_KIND, time, ierr))
370:   PetscCallA(DMSetOutputSequenceNumber(dmA, 0_PETSC_INT_KIND, time, ierr))

372:   PetscCallA(VecView(U, viewer, ierr))
373:   PetscCallA(VecView(A, viewer, ierr))

375:   ! Saving U and Alpha in one shot.
376:   ! For this, we need to cheat and change the Vec's name
377:   ! Note that in the end we write variables one component at a time,
378:   ! so that there is no real value in doing this
379:   PetscCallA(DMSetOutputSequenceNumber(dmUA, 1_PETSC_INT_KIND, time, ierr))
380:   PetscCallA(DMGetGlobalVector(dmUA, tmpVec, ierr))
381:   PetscCallA(VecCopy(UA, tmpVec, ierr))
382:   PetscCallA(PetscObjectSetName(tmpVec, 'U', ierr))
383:   PetscCallA(VecView(tmpVec, viewer, ierr))

385:   ! Reading nodal variables in Exodus file
386:   PetscCallA(VecSet(tmpVec, -1000.0_PETSC_REAL_KIND, ierr))
387:   PetscCallA(VecLoad(tmpVec, viewer, ierr))
388:   PetscCallA(VecAXPY(UA, -1.0_PETSC_REAL_KIND, tmpVec, ierr))
389:   PetscCallA(VecNorm(UA, NORM_INFINITY, norm, ierr))
390:   if (norm > PETSC_SQRT_MACHINE_EPSILON) then
391:     write (IOBuffer, '("UAlpha ||Vin - Vout|| = ",ES12.5)') norm
392:   end if
393:   PetscCallA(DMRestoreGlobalVector(dmUA, tmpVec, ierr))

395:   ! same thing with the UA2 Vec obtained from the superDM
396:   PetscCallA(DMGetGlobalVector(dmUA2, tmpVec, ierr))
397:   PetscCallA(VecCopy(UA2, tmpVec, ierr))
398:   PetscCallA(PetscObjectSetName(tmpVec, 'U', ierr))
399:   PetscCallA(DMSetOutputSequenceNumber(dmUA2, 2_PETSC_INT_KIND, time, ierr))
400:   PetscCallA(VecView(tmpVec, viewer, ierr))

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

412:   ! Building and saving Sigma
413:   !   We set sigma_0 = rank (to see partitioning)
414:   !          sigma_1 = cell set ID
415:   !          sigma_2 = x_coordinate of the cell center of mass
416:   PetscCallA(DMGetCoordinateSection(dmS, coordSection, ierr))
417:   PetscCallA(DMGetCoordinatesLocal(dmS, coord, ierr))
418:   PetscCallA(DMGetLabelIdIS(dmS, 'Cell Sets', csIS, ierr))
419:   PetscCallA(DMGetLabelSize(dmS, 'Cell Sets', numCS, ierr))
420:   PetscCallA(ISGetIndices(csIS, csID, ierr))

422:   do set = 1, numCS
423:     PetscCallA(DMGetStratumIS(dmS, 'Cell Sets', csID(set), cellIS, ierr))
424:     PetscCallA(ISGetIndices(cellIS, cellID, ierr))
425:     PetscCallA(ISGetSize(cellIS, numCells, ierr))
426:     do cell = 1, numCells
427:       PetscCallA(DMPlexVecGetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), PETSC_NULL_INTEGER, cval, ierr))
428:       PetscCallA(DMPlexVecGetClosure(dmS, coordSection, coord, cellID(cell), PETSC_NULL_INTEGER, xyz, ierr))
429:       cval(1) = rank
430:       cval(2) = csID(set)
431:       cval(3) = 0.0_PETSC_REAL_KIND
432:       do c = 1, size(xyz), sdim
433:         cval(3) = cval(3) + xyz(c)
434:       end do
435:       cval(3) = cval(3)*sdim/size(xyz)
436:       PetscCallA(DMPlexVecSetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval, INSERT_ALL_VALUES, ierr))
437:       PetscCallA(DMPlexVecRestoreClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), PETSC_NULL_INTEGER, cval, ierr))
438:       PetscCallA(DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID(cell), PETSC_NULL_INTEGER, xyz, ierr))
439:     end do
440:     PetscCallA(ISRestoreIndices(cellIS, cellID, ierr))
441:     PetscCallA(ISDestroy(cellIS, ierr))
442:   end do
443:   PetscCallA(ISRestoreIndices(csIS, csID, ierr))
444:   PetscCallA(ISDestroy(csIS, ierr))
445:   PetscCallA(VecViewFromOptions(S, PETSC_NULL_OBJECT, '-s_vec_view', ierr))

447:   ! Writing zonal variables in Exodus file
448:   PetscCallA(DMSetOutputSequenceNumber(dmS, 0_PETSC_INT_KIND, time, ierr))
449:   PetscCallA(VecView(S, viewer, ierr))

451:   ! Reading zonal variables in Exodus file */
452:   PetscCallA(DMGetGlobalVector(dmS, tmpVec, ierr))
453:   PetscCallA(VecSet(tmpVec, -1000.0_PETSC_REAL_KIND, ierr))
454:   PetscCallA(PetscObjectSetName(tmpVec, 'Sigma', ierr))
455:   PetscCallA(VecLoad(tmpVec, viewer, ierr))
456:   PetscCallA(VecAXPY(S, -1.0_PETSC_REAL_KIND, tmpVec, ierr))
457:   PetscCallA(VecNorm(S, NORM_INFINITY, norm, ierr))
458:   if (norm > PETSC_SQRT_MACHINE_EPSILON) then
459:     write (IOBuffer, '("Sigma ||Vin - Vout|| = ",ES12.5)') norm
460:   end if
461:   PetscCallA(DMRestoreGlobalVector(dmS, tmpVec, ierr))

463:   PetscCallA(DMRestoreGlobalVector(dmUA2, UA2, ierr))
464:   PetscCallA(DMRestoreGlobalVector(dmUA, UA, ierr))
465:   PetscCallA(DMRestoreGlobalVector(dmS, S, ierr))
466:   PetscCallA(DMRestoreGlobalVector(dmA, A, ierr))
467:   PetscCallA(DMRestoreGlobalVector(dmU, U, ierr))
468:   PetscCallA(DMRestoreGlobalVector(pdm, X, ierr))
469:   PetscCallA(DMDestroy(dmU, ierr))
470:   PetscCallA(ISDestroy(isU, ierr))
471:   PetscCallA(DMDestroy(dmA, ierr))
472:   PetscCallA(ISDestroy(isA, ierr))
473:   PetscCallA(DMDestroy(dmS, ierr))
474:   PetscCallA(ISDestroy(isS, ierr))
475:   PetscCallA(DMDestroy(dmUA, ierr))
476:   PetscCallA(ISDestroy(isUA, ierr))
477:   PetscCallA(DMDestroy(dmUA2, ierr))
478:   PetscCallA(DMDestroy(pdm, ierr))
479:   if (numProc > 1) then
480:     PetscCallA(DMDestroy(dm, ierr))
481:   end if

483:   deallocate (pStartDepth)
484:   deallocate (pEndDepth)

486:   PetscCallA(PetscViewerDestroy(viewer, ierr))
487:   PetscCallA(PetscFinalize(ierr))
488: end program ex62f90

490: ! /*TEST
491: !
492: ! build:
493: !   requires: exodusii pnetcdf !complex
494: !   # 2D seq
495: ! test:
496: !   suffix: 0
497: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 1
498: !   #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
499: ! test:
500: !   suffix: 1
501: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 1
502: !
503: ! test:
504: !   suffix: 2
505: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 1
506: !   #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
507: ! test:
508: !   suffix: 3
509: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 2
510: ! test:
511: !   suffix: 4
512: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 2
513: ! test:
514: !   suffix: 5
515: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 2
516: !   # 2D par
517: ! test:
518: !   suffix: 6
519: !   nsize: 2
520: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_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: 7
524: !   nsize: 2
525: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 1
526: ! test:
527: !   suffix: 8
528: !   nsize: 2
529: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 1
530: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name
531: ! test:
532: !   suffix: 9
533: !   nsize: 2
534: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -o FourSquareT-large_out.exo -dm_view -petscpartitioner_type simple -order 2
535: ! test:
536: !   suffix: 10
537: !   nsize: 2
538: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -o FourSquareQ-large_out.exo -dm_view -petscpartitioner_type simple -order 2
539: ! test:
540: !   # Something is now broken with parallel read/write for whatever shape H is
541: !   TODO: broken
542: !   suffix: 11
543: !   nsize: 2
544: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -o FourSquareH-large_out.exo -dm_view -petscpartitioner_type simple -order 2
545: !   #3d seq
546: ! test:
547: !   suffix: 12
548: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 1
549: ! test:
550: !   suffix: 13
551: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 1
552: !   #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
553: ! test:
554: !   suffix: 14
555: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 2
556: ! test:
557: !   suffix: 15
558: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 2
559: !   #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
560: !   #3d par
561: ! test:
562: !   suffix: 16
563: !   nsize: 2
564: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 1
565: ! test:
566: !   suffix: 17
567: !   nsize: 2
568: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 1
569: !   #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
570: ! test:
571: !   suffix: 18
572: !   nsize: 2
573: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -o FourBrickHex-large_out.exo -dm_view -petscpartitioner_type simple -order 2
574: ! test:
575: !   suffix: 19
576: !   nsize: 2
577: !   args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -o FourBrickTet-large_out.exo -dm_view -petscpartitioner_type simple -order 2
578: !   #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
579: !
580: ! TEST*/