Actual source code: ex26f90.F90

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

  7:     ! Get the Fortran kind associated with PetscInt and PetscReal so that we can use literal constants.
  8:     PetscInt                           :: dummyPetscInt
  9:     PetscReal                          :: dummyPetscreal
 10:     integer,parameter                  :: kPI = kind(dummyPetscInt)
 11:     integer,parameter                  :: kPR = kind(dummyPetscReal)

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

 36:     Character(len=MXSTLN)              :: sJunk
 37:     PetscInt                           :: numstep = 3, step
 38:     Integer                            :: numNodalVar,numZonalVar
 39:     character(len=MXNAME),dimension(4) :: nodalVarName2D = ["U_x  ", &
 40:                                                             "U_y  ", &
 41:                                                             "Alpha", &
 42:                                                             "Beta "]
 43:     character(len=MXNAME),dimension(3) :: zonalVarName2D = ["Sigma_11", &
 44:                                                             "Sigma_12", &
 45:                                                             "Sigma_22"]
 46:     character(len=MXNAME),dimension(5) :: nodalVarName3D = ["U_x  ", &
 47:                                                             "U_y  ", &
 48:                                                             "U_z  ", &
 49:                                                             "Alpha", &
 50:                                                             "Beta "]
 51:     character(len=MXNAME),dimension(6) :: zonalVarName3D = ["Sigma_11", &
 52:                                                             "Sigma_22", &
 53:                                                             "Sigma_33", &
 54:                                                             "Sigma_23", &
 55:                                                             "Sigma_13", &
 56:                                                             "Sigma_12"]
 57:     logical,dimension(:,:),pointer     :: truthtable

 59:     type(tIS)                          :: cellIS
 60:     PetscInt,dimension(:),pointer      :: cellID
 61:     PetscInt                           :: numCells, cell, closureSize
 62:     PetscInt,dimension(:),pointer      :: closureA,closure

 64:     type(tPetscSection)                :: sectionUA,coordSection
 65:     type(tVec)                         :: UALoc,coord
 66:     PetscScalar,dimension(:),pointer   :: cval,xyz
 67:     PetscInt                           :: dofUA,offUA,c

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

 90:     type(tPetscSF)                      :: migrationSF
 91:     PetscPartitioner                    :: part
 92:     PetscLayout                         :: layout

 94:     type(tVec)                          :: tmpVec
 95:     PetscReal                           :: norm
 96:     PetscReal                           :: time = 1.234_kPR

 98:     PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER,ierr))
 99:     if (ierr /= 0) then
100:       print*,'Unable to initialize PETSc'
101:       stop
102:     endif

104:     PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
105:     PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,numProc,ierr))
106:     PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-i',ifilename,flg,ierr))
107:     PetscCheckA(flg,PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,'missing input file name -i ')
108:     PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-o',ofilename,flg,ierr))
109:     PetscCheckA(flg,PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,'missing output file name -o ')
110:     PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-order',order,flg,ierr))
111:     if ((order > 2) .or. (order < 1)) then
112:         write(IOBuffer,'("Unsupported polynomial order ", I2, " not in [1,2]")') order
113:         SETERRA(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,IOBuffer)
114:     end if

116:     ! Read the mesh in any supported format
117:     PetscCallA(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename,PETSC_NULL_CHARACTER,PETSC_TRUE,dm,ierr))
118:     PetscCallA(DMPlexDistributeSetDefault(dm,PETSC_FALSE,ierr));
119:     PetscCallA(DMSetFromOptions(dm,ierr));
120:     PetscCallA(DMGetDimension(dm, sdim,ierr))
121:     PetscCallA(DMViewFromOptions(dm, PETSC_NULL_OPTIONS,'-dm_view',ierr));

123:     ! Create the exodus result file

125:     ! enable exodus debugging information
126:     PetscCallA(exopts(EXVRBS+EXDEBG,ierr))
127:     ! Create the exodus file
128:     PetscCallA(PetscViewerExodusIIOpen(PETSC_COMM_WORLD,ofilename,FILE_MODE_WRITE,viewer,ierr))
129:     ! The long way would be
130:     !
131:     ! PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD,viewer,ierr))
132:     ! PetscCallA(PetscViewerSetType(viewer,PETSCVIEWEREXODUSII,ierr))
133:     ! PetscCallA(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE,ierr))
134:     ! PetscCallA(PetscViewerFileSetName(viewer,ofilename,ierr))

136:     ! set the mesh order
137:     PetscCallA(PetscViewerExodusIISetOrder(viewer,order,ierr))
138:     PetscCallA(PetscViewerView(viewer,PETSC_VIEWER_STDOUT_WORLD,ierr))
139:     !
140:     !    Notice how the exodus file is actually NOT open at this point (exoid is -1)
141:     !    Since we are overwriting the file (mode is FILE_MODE_WRITE), we are going to have to
142:     !    write the geometry (the DM), which can only be done on a brand new file.
143:     !

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

180:     !    An exodusII truth table specifies which fields are saved at which time step
181:     !    It speeds up I/O but reserving space for fields in the file ahead of time.
182:     PetscCallA(PetscViewerExodusIIGetId(viewer,exoid,ierr))
183:     PetscCallA(exinq(exoid, EX_INQ_ELEM_BLK,numCS,PETSC_NULL_REAL,sjunk,ierr))
184:     allocate(truthtable(numCS,numZonalVar))
185:     truthtable = .true.
186:     PetscCallA(expvtt(exoid, numCS, numZonalVar, truthtable, ierr))
187:     deallocate(truthtable)

189:     !   Writing time step information in the file. Note that this is currently broken in the exodus library for netcdf4 (HDF5-based) files */
190:     do step = 1,numstep
191:         PetscCallA(exptim(exoid,step,Real(step,kind=kPR),ierr))
192:     end do

194:     PetscCallA(PetscObjectGetComm(dm,comm,ierr))
195:     PetscCallA(PetscSectionCreate(comm, section,ierr))
196:     PetscCallA(PetscSectionSetNumFields(section, 3_kPI,ierr))
197:     PetscCallA(PetscSectionSetFieldName(section, fieldU, 'U',ierr))
198:     PetscCallA(PetscSectionSetFieldName(section, fieldA, 'Alpha',ierr))
199:     PetscCallA(PetscSectionSetFieldName(section, fieldS, 'Sigma',ierr))
200:     PetscCallA(DMPlexGetChart(dm, pStart, pEnd,ierr))
201:     PetscCallA(PetscSectionSetChart(section, pStart, pEnd,ierr))

203:     allocate(pStartDepth(sdim+1))
204:     allocate(pEndDepth(sdim+1))
205:     do d = 1, sdim+1
206:         PetscCallA(DMPlexGetDepthStratum(dm, d-1, pStartDepth(d), pEndDepth(d),ierr))
207:     end do

209:     ! Vector field U, Scalar field Alpha, Tensor field Sigma
210:     PetscCallA(PetscSectionSetFieldComponents(section, fieldU, sdim,ierr));
211:     PetscCallA(PetscSectionSetFieldComponents(section, fieldA, 1_kPI,ierr));
212:     PetscCallA(PetscSectionSetFieldComponents(section, fieldS, sdim*(sdim+1)/2,ierr));

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

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

308:     PetscCallA(DMSetUseNatural(dm,PETSC_TRUE,ierr))
309:     PetscCallA(DMPlexGetPartitioner(dm,part,ierr))
310:     PetscCallA(PetscPartitionerSetFromOptions(part,ierr))
311:     PetscCallA(DMPlexDistribute(dm,0_kPI,migrationSF,pdm,ierr))

313:     if (numProc > 1) then
314:         PetscCallA(DMPlexSetMigrationSF(pdm,migrationSF,ierr))
315:         PetscCallA(PetscSFDestroy(migrationSF,ierr))
316:     else
317:         pdm = dm
318:     end if
319:     PetscCallA(DMViewFromOptions(pdm,PETSC_NULL_OPTIONS,'-dm_view',ierr))

321:     ! Get DM and IS for each field of dm
322:     PetscCallA(DMCreateSubDM(pdm, 1_kPI, fieldU,  isU,  dmU,ierr))
323:     PetscCallA(DMCreateSubDM(pdm, 1_kPI, fieldA,  isA,  dmA,ierr))
324:     PetscCallA(DMCreateSubDM(pdm, 1_kPI, fieldS,  isS,  dmS,ierr))
325:     PetscCallA(DMCreateSubDM(pdm, 2_kPI, fieldUA, isUA, dmUA,ierr))

327:     !Create the exodus result file
328:     allocate(dmList(2))
329:     dmList(1) = dmU;
330:     dmList(2) = dmA;
331:     PetscCallA(DMCreateSuperDM(dmList,2_kPI,PETSC_NULL_IS,dmUA2,ierr))
332:     deallocate(dmList)

334:     PetscCallA(DMGetGlobalVector(pdm,   X,ierr))
335:     PetscCallA(DMGetGlobalVector(dmU,  U,ierr))
336:     PetscCallA(DMGetGlobalVector(dmA,  A,ierr))
337:     PetscCallA(DMGetGlobalVector(dmS,  S,ierr))
338:     PetscCallA(DMGetGlobalVector(dmUA, UA,ierr))
339:     PetscCallA(DMGetGlobalVector(dmUA2, UA2,ierr))

341:     PetscCallA(PetscObjectSetName(U,  'U',ierr))
342:     PetscCallA(PetscObjectSetName(A,  'Alpha',ierr))
343:     PetscCallA(PetscObjectSetName(S,  'Sigma',ierr))
344:     PetscCallA(PetscObjectSetName(UA, 'UAlpha',ierr))
345:     PetscCallA(PetscObjectSetName(UA2, 'UAlpha2',ierr))
346:     PetscCallA(VecSet(X, -111.0_kPR,ierr))

348:     ! 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 */
349:     PetscCallA(DMGetLocalSection(dmUA, sectionUA,ierr))
350:     PetscCallA(DMGetLocalVector(dmUA, UALoc,ierr))
351:     PetscCallA(VecGetArrayF90(UALoc, cval,ierr))
352:     PetscCallA(DMGetCoordinateSection(dmUA, coordSection,ierr))
353:     PetscCallA(DMGetCoordinatesLocal(dmUA, coord,ierr))
354:     PetscCallA(DMPlexGetChart(dmUA, pStart, pEnd,ierr))

356:     do p = pStart,pEnd-1
357:         PetscCallA(PetscSectionGetDof(sectionUA, p, dofUA,ierr))
358:         if (dofUA > 0) then
359:             PetscCallA(PetscSectionGetOffset(sectionUA, p, offUA,ierr))
360:             PetscCallA(DMPlexVecGetClosure(dmUA, coordSection, coord, p, xyz,ierr))
361:             closureSize = size(xyz)
362:             do i = 1,sdim
363:                 do j = 0, closureSize-1,sdim
364:                     cval(offUA+i) = cval(offUA+i) + xyz(j/sdim+i)
365:                 end do
366:                 cval(offUA+i) = cval(offUA+i) * sdim / closureSize;
367:                 cval(offUA+sdim+1) = cval(offUA+sdim+1) + cval(offUA+i)**2
368:             end do
369:             PetscCallA(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, xyz,ierr))
370:         end if
371:     end do

373:     PetscCallA(VecRestoreArrayF90(UALoc, cval,ierr))
374:     PetscCallA(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA,ierr))
375:     PetscCallA(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA,ierr))
376:     PetscCallA(DMRestoreLocalVector(dmUA, UALoc,ierr))

378:     !Update X
379:     PetscCallA(VecISCopy(X, isUA, SCATTER_FORWARD, UA,ierr))
380:     ! Restrict to U and Alpha
381:     PetscCallA(VecISCopy(X, isU, SCATTER_REVERSE, U,ierr))
382:     PetscCallA(VecISCopy(X, isA, SCATTER_REVERSE, A,ierr))
383:     PetscCallA(VecViewFromOptions(UA, PETSC_NULL_OPTIONS, '-ua_vec_view',ierr))
384:     PetscCallA(VecViewFromOptions(U, PETSC_NULL_OPTIONS, '-u_vec_view',ierr))
385:     PetscCallA(VecViewFromOptions(A, PETSC_NULL_OPTIONS, '-a_vec_view',ierr))
386:     ! restrict to UA2
387:     PetscCallA(VecISCopy(X, isUA, SCATTER_REVERSE, UA2,ierr))
388:     PetscCallA(VecViewFromOptions(UA2, PETSC_NULL_OPTIONS, '-ua2_vec_view',ierr))

390:     ! Getting Natural Vec
391:     PetscCallA(DMSetOutputSequenceNumber(dmU, 0_kPI, time, ierr))
392:     PetscCallA(DMSetOutputSequenceNumber(dmA, 0_kPI, time, ierr))

394:     PetscCallA(VecView(U, viewer,ierr))
395:     PetscCallA(VecView(A, viewer,ierr))

397:     !Saving U and Alpha in one shot.
398:     !For this, we need to cheat and change the Vec's name
399:     !Note that in the end we write variables one component at a time,
400:     !so that there is no real value in doing this
401:     PetscCallA(DMSetOutputSequenceNumber(dmUA,1_kPI,time,ierr))
402:     PetscCallA(DMGetGlobalVector(dmUA, tmpVec,ierr))
403:     PetscCallA(VecCopy(UA, tmpVec,ierr))
404:     PetscCallA(PetscObjectSetName(tmpVec, 'U',ierr))
405:     PetscCallA(VecView(tmpVec, viewer,ierr))

407:     ! Reading nodal variables in Exodus file
408:     PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr))
409:     PetscCallA(VecLoad(tmpVec, viewer,ierr))
410:     PetscCallA(VecAXPY(UA, -1.0_kPR, tmpVec,ierr))
411:     PetscCallA(VecNorm(UA, NORM_INFINITY, norm,ierr))
412:     if (norm > PETSC_SQRT_MACHINE_EPSILON) then
413:         write(IOBuffer,'("UAlpha ||Vin - Vout|| = ",ES12.5)') norm
414:     end if
415:     PetscCallA(DMRestoreGlobalVector(dmUA, tmpVec,ierr))

417:     ! same thing with the UA2 Vec obtained from the superDM
418:     PetscCallA(DMGetGlobalVector(dmUA2, tmpVec,ierr))
419:     PetscCallA(VecCopy(UA2, tmpVec,ierr))
420:     PetscCallA(PetscObjectSetName(tmpVec, 'U',ierr))
421:     PetscCallA(DMSetOutputSequenceNumber(dmUA2,2_kPI,time,ierr))
422:     PetscCallA(VecView(tmpVec, viewer,ierr))

424:     ! Reading nodal variables in Exodus file
425:     PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr))
426:     PetscCallA(VecLoad(tmpVec,viewer,ierr))
427:     PetscCallA(VecAXPY(UA2, -1.0_kPR, tmpVec,ierr))
428:     PetscCallA(VecNorm(UA2, NORM_INFINITY, norm,ierr))
429:     if (norm > PETSC_SQRT_MACHINE_EPSILON) then
430:         write(IOBuffer,'("UAlpha2 ||Vin - Vout|| = ",ES12.5)') norm
431:     end if
432:     PetscCallA(DMRestoreGlobalVector(dmUA2, tmpVec,ierr))

434:     ! Building and saving Sigma
435:     !   We set sigma_0 = rank (to see partitioning)
436:     !          sigma_1 = cell set ID
437:     !          sigma_2 = x_coordinate of the cell center of mass
438:     PetscCallA(DMGetCoordinateSection(dmS, coordSection,ierr))
439:     PetscCallA(DMGetCoordinatesLocal(dmS, coord,ierr))
440:     PetscCallA(DMGetLabelIdIS(dmS, 'Cell Sets', csIS,ierr))
441:     PetscCallA(DMGetLabelSize(dmS, 'Cell Sets',numCS,ierr))
442:     PetscCallA(ISGetIndicesF90(csIS, csID,ierr))

444:     do set = 1, numCS
445:         PetscCallA(DMGetStratumIS(dmS, 'Cell Sets', csID(set), cellIS,ierr))
446:         PetscCallA(ISGetIndicesF90(cellIS, cellID,ierr))
447:         PetscCallA(ISGetSize(cellIS, numCells,ierr))
448:         do cell = 1,numCells
449:             PetscCallA(DMPlexVecGetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval,ierr))
450:             PetscCallA(DMPlexVecGetClosure(dmS, coordSection, coord, cellID(cell), xyz,ierr))
451:             cval(1) = rank
452:             cval(2) = csID(set)
453:             cval(3) = 0.0_kPR
454:             do c = 1, size(xyz),sdim
455:                 cval(3) = cval(3) + xyz(c)
456:             end do
457:             cval(3) = cval(3) * sdim / size(xyz)
458:             PetscCallA(DMPlexVecSetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval, INSERT_ALL_VALUES,ierr))
459:             PetscCallA(DMPlexVecRestoreClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval,ierr))
460:             PetscCallA(DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID(cell), xyz,ierr))
461:         end do
462:         PetscCallA(ISRestoreIndicesF90(cellIS, cellID,ierr))
463:         PetscCallA(ISDestroy(cellIS,ierr))
464:     end do
465:     PetscCallA(ISRestoreIndicesF90(csIS, csID,ierr))
466:     PetscCallA(ISDestroy(csIS,ierr))
467:     PetscCallA(VecViewFromOptions(S, PETSC_NULL_OPTIONS, '-s_vec_view',ierr))

469:     ! Writing zonal variables in Exodus file
470:     PetscCallA(DMSetOutputSequenceNumber(dmS,0_kPI,time,ierr))
471:     PetscCallA(VecView(S,viewer,ierr))

473:     ! Reading zonal variables in Exodus file */
474:     PetscCallA(DMGetGlobalVector(dmS, tmpVec,ierr))
475:     PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr))
476:     PetscCallA(PetscObjectSetName(tmpVec, 'Sigma',ierr))
477:     PetscCallA(VecLoad(tmpVec,viewer,ierr))
478:     PetscCallA(VecAXPY(S, -1.0_kPR, tmpVec,ierr))
479:     PetscCallA(VecNorm(S, NORM_INFINITY,norm,ierr))
480:     if (norm > PETSC_SQRT_MACHINE_EPSILON) then
481:        write(IOBuffer,'("Sigma ||Vin - Vout|| = ",ES12.5)') norm
482:     end if
483:     PetscCallA(DMRestoreGlobalVector(dmS, tmpVec,ierr))

485:     PetscCallA(DMRestoreGlobalVector(dmUA2, UA2,ierr))
486:     PetscCallA(DMRestoreGlobalVector(dmUA, UA,ierr))
487:     PetscCallA(DMRestoreGlobalVector(dmS,  S,ierr))
488:     PetscCallA(DMRestoreGlobalVector(dmA,  A,ierr))
489:     PetscCallA(DMRestoreGlobalVector(dmU,  U,ierr))
490:     PetscCallA(DMRestoreGlobalVector(pdm,   X,ierr))
491:     PetscCallA(DMDestroy(dmU,ierr))
492:     PetscCallA(ISDestroy(isU,ierr))
493:     PetscCallA(DMDestroy(dmA,ierr))
494:     PetscCallA(ISDestroy(isA,ierr))
495:     PetscCallA(DMDestroy(dmS,ierr))
496:     PetscCallA(ISDestroy(isS,ierr))
497:     PetscCallA(DMDestroy(dmUA,ierr))
498:     PetscCallA(ISDestroy(isUA,ierr))
499:     PetscCallA(DMDestroy(dmUA2,ierr))
500:     PetscCallA(DMDestroy(pdm,ierr))
501:     if (numProc > 1) then
502:         PetscCallA(DMDestroy(dm,ierr))
503:     end if

505:     deallocate(pStartDepth)
506:     deallocate(pEndDepth)

508:     PetscCallA(PetscViewerDestroy(viewer,ierr))
509:     PetscCallA(PetscFinalize(ierr))
510: end program ex62f90

512: ! /*TEST
513: !
514: ! build:
515: !   requires: exodusii pnetcdf !complex
516: !   # 2D seq
517: ! test:
518: !   suffix: 0
519: !   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
520: !   #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
521: ! test:
522: !   suffix: 1
523: !   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
524: !
525: ! test:
526: !   suffix: 2
527: !   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
528: !   #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
529: ! test:
530: !   suffix: 3
531: !   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
532: ! test:
533: !   suffix: 4
534: !   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
535: ! test:
536: !   suffix: 5
537: !   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
538: !   # 2D par
539: ! test:
540: !   suffix: 6
541: !   nsize: 2
542: !   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
543: !   #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
544: ! test:
545: !   suffix: 7
546: !   nsize: 2
547: !   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
548: ! test:
549: !   suffix: 8
550: !   nsize: 2
551: !   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
552: !   #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name
553: ! test:
554: !   suffix: 9
555: !   nsize: 2
556: !   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
557: ! test:
558: !   suffix: 10
559: !   nsize: 2
560: !   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
561: ! test:
562: !   # Something is now broken with parallel read/write for whatever shape H is
563: !   TODO: broken
564: !   suffix: 11
565: !   nsize: 2
566: !   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

568: !   #3d seq
569: ! test:
570: !   suffix: 12
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 1
572: ! test:
573: !   suffix: 13
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 1
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: ! test:
577: !   suffix: 14
578: !   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
579: ! test:
580: !   suffix: 15
581: !   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
582: !   #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
583: !   #3d par
584: ! test:
585: !   suffix: 16
586: !   nsize: 2
587: !   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
588: ! test:
589: !   suffix: 17
590: !   nsize: 2
591: !   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
592: !   #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
593: ! test:
594: !   suffix: 18
595: !   nsize: 2
596: !   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
597: ! test:
598: !   suffix: 19
599: !   nsize: 2
600: !   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
601: !   #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

603: ! TEST*/