1program ex62f90 2#include "petsc/finclude/petsc.h" 3 use petsc 4 implicit none 5#include "exodusII.inc" 6 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) 12 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 PetscInt :: sdim,d,pStart,pEnd,p,numCS,set,i,j 29 PetscMPIInt :: rank,numProc 30 PetscBool :: flg 31 PetscErrorCode :: ierr 32 MPI_Comm :: comm 33 type(tPetscViewer) :: viewer 34 35 Character(len=MXSTLN) :: sJunk 36 PetscInt :: numstep = 3, step 37 PetscInt :: numNodalVar,numZonalVar 38 character(len=MXSTLN) :: nodalVarName(4) 39 character(len=MXSTLN) :: zonalVarName(6) 40 logical,dimension(:,:),pointer :: truthtable 41 42 type(tIS) :: cellIS 43 PetscInt,dimension(:),pointer :: cellID 44 PetscInt :: numCells, cell, closureSize 45 PetscInt,dimension(:),pointer :: closureA,closure 46 47 type(tPetscSection) :: sectionUA,coordSection 48 type(tVec) :: UALoc,coord 49 PetscScalar,dimension(:),pointer :: cval,xyz 50 PetscInt :: dofUA,offUA,c 51 52 ! dof layout ordered by increasing height in the DAG: cell, face, edge, vertex 53 PetscInt,dimension(3),target :: dofS2D = [0, 0, 3] 54 PetscInt,dimension(3),target :: dofUP1Tri = [2, 0, 0] 55 PetscInt,dimension(3),target :: dofAP1Tri = [1, 0, 0] 56 PetscInt,dimension(3),target :: dofUP2Tri = [2, 2, 0] 57 PetscInt,dimension(3),target :: dofAP2Tri = [1, 1, 0] 58 PetscInt,dimension(3),target :: dofUP1Quad = [2, 0, 0] 59 PetscInt,dimension(3),target :: dofAP1Quad = [1, 0, 0] 60 PetscInt,dimension(3),target :: dofUP2Quad = [2, 2, 2] 61 PetscInt,dimension(3),target :: dofAP2Quad = [1, 1, 1] 62 PetscInt,dimension(4),target :: dofS3D = [0, 0, 0, 6] 63 PetscInt,dimension(4),target :: dofUP1Tet = [3, 0, 0, 0] 64 PetscInt,dimension(4),target :: dofAP1Tet = [1, 0, 0, 0] 65 PetscInt,dimension(4),target :: dofUP2Tet = [3, 3, 0, 0] 66 PetscInt,dimension(4),target :: dofAP2Tet = [1, 1, 0, 0] 67 PetscInt,dimension(4),target :: dofUP1Hex = [3, 0, 0, 0] 68 PetscInt,dimension(4),target :: dofAP1Hex = [1, 0, 0, 0] 69 PetscInt,dimension(4),target :: dofUP2Hex = [3, 3, 3, 3] 70 PetscInt,dimension(4),target :: dofAP2Hex = [1, 1, 1, 1] 71 PetscInt,dimension(:),pointer :: dofU,dofA,dofS 72 73 type(tPetscSF) :: migrationSF 74 PetscPartitioner :: part 75 76 type(tVec) :: tmpVec 77 PetscReal :: norm 78 PetscReal :: time = 1.234_kPR 79 80 PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER,ierr)) 81 if (ierr /= 0) then 82 print*,'Unable to initialize PETSc' 83 stop 84 endif 85 86 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)) 87 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,numProc,ierr)) 88 PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-i',ifilename,flg,ierr)) 89 PetscCheckA(flg,PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,'missing input file name -i <input file name>') 90 PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-o',ofilename,flg,ierr)) 91 PetscCheckA(flg,PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,'missing output file name -o <output file name>') 92 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-order',order,flg,ierr)) 93 if ((order > 2) .or. (order < 1)) then 94 write(IOBuffer,'("Unsupported polynomial order ", I2, " not in [1,2]")') order 95 SETERRA(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,IOBuffer) 96 end if 97 98 ! Read the mesh in any supported format 99 PetscCallA(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename,PETSC_NULL_CHARACTER,PETSC_TRUE,dm,ierr)) 100 PetscCallA(DMPlexDistributeSetDefault(dm,PETSC_FALSE,ierr)); 101 PetscCallA(DMSetFromOptions(dm,ierr)); 102 PetscCallA(DMGetDimension(dm, sdim,ierr)) 103 PetscCallA(DMViewFromOptions(dm, PETSC_NULL_OPTIONS,'-dm_view',ierr)); 104 105 ! Create the exodus result file 106 107 ! enable exodus debugging information 108 PetscCallA(exopts(EXVRBS+EXDEBG,ierr)) 109 ! Create the exodus file 110 PetscCallA(PetscViewerExodusIIOpen(PETSC_COMM_WORLD,ofilename,FILE_MODE_WRITE,viewer,ierr)) 111 ! The long way would be 112 ! 113 ! PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD,viewer,ierr)) 114 ! PetscCallA(PetscViewerSetType(viewer,PETSCVIEWEREXODUSII,ierr)) 115 ! PetscCallA(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE,ierr)) 116 ! PetscCallA(PetscViewerFileSetName(viewer,ofilename,ierr)) 117 118 ! set the mesh order 119 PetscCallA(PetscViewerExodusIISetOrder(viewer,order,ierr)) 120 PetscCallA(PetscViewerView(viewer,PETSC_VIEWER_STDOUT_WORLD,ierr)) 121 ! 122 ! Notice how the exodus file is actually NOT open at this point (exoid is -1) 123 ! Since we are overwriting the file (mode is FILE_MODE_WRITE), we are going to have to 124 ! write the geometry (the DM), which can only be done on a brand new file. 125 ! 126 127 ! Save the geometry to the file, erasing all previous content 128 PetscCallA(DMView(dm,viewer,ierr)) 129 PetscCallA(PetscViewerView(viewer,PETSC_VIEWER_STDOUT_WORLD,ierr)) 130 ! 131 ! Note how the exodus file is now open 132 ! 133 ! 'Format' the exodus result file, i.e. allocate space for nodal and zonal variables 134 select case(sdim) 135 case(2) 136 numNodalVar = 3 137 nodalVarName(1:numNodalVar) = ['U_x ','U_y ','Alpha'] 138 numZonalVar = 3 139 zonalVarName(1:numZonalVar) = ['Sigma_11','Sigma_22','Sigma_12'] 140 case(3) 141 numNodalVar = 4 142 nodalVarName(1:numNodalVar) = ['U_x ','U_y ','U_z ','Alpha'] 143 numZonalVar = 6 144 zonalVarName(1:numZonalVar) = ['Sigma_11','Sigma_22','Sigma_33','Sigma_23','Sigma_13','Sigma_12'] 145 case default 146 write(IOBuffer,'("No layout for dimension ",I2)') sdim 147 end select 148 PetscCallA(PetscViewerExodusIIGetId(viewer,exoid,ierr)) 149 PetscCallA(expvp(exoid, 'E', numZonalVar,ierr)) 150 PetscCallA(expvan(exoid, 'E', numZonalVar, zonalVarName,ierr)) 151 PetscCallA(expvp(exoid, 'N', numNodalVar,ierr)) 152 PetscCallA(expvan(exoid, 'N', numNodalVar, nodalVarName,ierr)) 153 PetscCallA(exinq(exoid, EX_INQ_ELEM_BLK,numCS,PETSC_NULL_REAL,sjunk,ierr)) 154 155 ! An exodusII truth table specifies which fields are saved at which time step 156 ! It speeds up I/O but reserving space for fields in the file ahead of time. 157 allocate(truthtable(numCS,numZonalVar)) 158 truthtable = .true. 159 PetscCallA(expvtt(exoid, numCS, numZonalVar, truthtable, ierr)) 160 deallocate(truthtable) 161 162 ! Writing time step information in the file. Note that this is currently broken in the exodus library for netcdf4 (HDF5-based) files */ 163 do step = 1,numstep 164 PetscCallA(exptim(exoid,step,Real(step,kind=kPR),ierr)) 165 end do 166 167 PetscCallA(PetscObjectGetComm(dm,comm,ierr)) 168 PetscCallA(PetscSectionCreate(comm, section,ierr)) 169 PetscCallA(PetscSectionSetNumFields(section, 3_kPI,ierr)) 170 PetscCallA(PetscSectionSetFieldName(section, fieldU, 'U',ierr)) 171 PetscCallA(PetscSectionSetFieldName(section, fieldA, 'Alpha',ierr)) 172 PetscCallA(PetscSectionSetFieldName(section, fieldS, 'Sigma',ierr)) 173 PetscCallA(DMPlexGetChart(dm, pStart, pEnd,ierr)) 174 PetscCallA(PetscSectionSetChart(section, pStart, pEnd,ierr)) 175 176 allocate(pStartDepth(sdim+1)) 177 allocate(pEndDepth(sdim+1)) 178 do d = 1, sdim+1 179 PetscCallA(DMPlexGetDepthStratum(dm, d-1, pStartDepth(d), pEndDepth(d),ierr)) 180 end do 181 182 ! Vector field U, Scalar field Alpha, Tensor field Sigma 183 PetscCallA(PetscSectionSetFieldComponents(section, fieldU, sdim,ierr)); 184 PetscCallA(PetscSectionSetFieldComponents(section, fieldA, 1_kPI,ierr)); 185 PetscCallA(PetscSectionSetFieldComponents(section, fieldS, sdim*(sdim+1)/2,ierr)); 186 187 ! Going through cell sets then cells, and setting up storage for the sections 188 PetscCallA(DMGetLabelSize(dm, 'Cell Sets', numCS, ierr)) 189 PetscCallA(DMGetLabelIdIS(dm, 'Cell Sets', csIS, ierr)) 190 PetscCallA(ISGetIndicesF90(csIS, csID, ierr)) 191 do set = 1,numCS 192 !!! This is pointless but will avoid a warning with gfortran and -Werror=maybe-uninitialized 193 nullify(dofA) 194 nullify(dofU) 195 nullify(dofS) 196 PetscCallA(DMGetStratumSize(dm, 'Cell Sets', csID(set), numCells,ierr)) 197 PetscCallA(DMGetStratumIS(dm, 'Cell Sets', csID(set), cellIS,ierr)) 198 if (numCells > 0) then 199 select case(sdim) 200 case(2) 201 dofs => dofS2D 202 case(3) 203 dofs => dofS3D 204 case default 205 write(IOBuffer,'("No layout for dimension ",I2)') sdim 206 SETERRA(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,IOBuffer) 207 end select ! sdim 208 209 ! Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes 210 ! It will not be enough to identify more exotic elements like pyramid or prisms... */ 211 PetscCallA(ISGetIndicesF90(cellIS, cellID,ierr)) 212 nullify(closureA) 213 PetscCallA(DMPlexGetTransitiveClosure(dm,cellID(1), PETSC_TRUE, closureA,ierr)) 214 select case(size(closureA)/2) 215 case(7) ! Tri 216 if (order == 1) then 217 dofU => dofUP1Tri 218 dofA => dofAP1Tri 219 else 220 dofU => dofUP2Tri 221 dofA => dofAP2Tri 222 end if 223 case(9) ! Quad 224 if (order == 1) then 225 dofU => dofUP1Quad 226 dofA => dofAP1Quad 227 else 228 dofU => dofUP2Quad 229 dofA => dofAP2Quad 230 end if 231 case(15) ! Tet 232 if (order == 1) then 233 dofU => dofUP1Tet 234 dofA => dofAP1Tet 235 else 236 dofU => dofUP2Tet 237 dofA => dofAP2Tet 238 end if 239 case(27) ! Hex 240 if (order == 1) then 241 dofU => dofUP1Hex 242 dofA => dofAP1Hex 243 else 244 dofU => dofUP2Hex 245 dofA => dofAP2Hex 246 end if 247 case default 248 write(IOBuffer,'("Unknown element with closure size ",I2)') size(closureA)/2 249 SETERRA(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,IOBuffer) 250 end select 251 PetscCallA(DMPlexRestoreTransitiveClosure(dm, cellID(1), PETSC_TRUE,closureA,ierr)) 252 do cell = 1,numCells! 253 nullify(closure) 254 PetscCallA(DMPlexGetTransitiveClosure(dm, cellID(cell), PETSC_TRUE, closure,ierr)) 255 do p = 1,size(closure),2 256 ! find the depth of p 257 do d = 1,sdim+1 258 if ((closure(p) >= pStartDepth(d)) .and. (closure(p) < pEndDepth(d))) then 259 PetscCallA(PetscSectionSetDof(section, closure(p), dofU(d)+dofA(d)+dofS(d),ierr)) 260 PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldU, dofU(d),ierr)) 261 PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldA, dofA(d),ierr)) 262 PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldS, dofS(d),ierr)) 263 end if ! closure(p) 264 end do ! d 265 end do ! p 266 PetscCallA(DMPlexRestoreTransitiveClosure(dm, cellID(cell), PETSC_TRUE, closure,ierr)) 267 end do ! cell 268 PetscCallA(ISRestoreIndicesF90(cellIS, cellID,ierr)) 269 PetscCallA(ISDestroy(cellIS,ierr)) 270 end if ! numCells 271 end do ! set 272 PetscCallA(ISRestoreIndicesF90(csIS, csID,ierr)) 273 PetscCallA(ISDestroy(csIS,ierr)) 274 PetscCallA(PetscSectionSetUp(section,ierr)) 275 PetscCallA(DMSetLocalSection(dm, section,ierr)) 276 PetscCallA(PetscObjectViewFromOptions(section, PETSC_NULL_SECTION, '-dm_section_view',ierr)) 277 PetscCallA(PetscSectionDestroy(section,ierr)) 278 279 PetscCallA(DMSetUseNatural(dm,PETSC_TRUE,ierr)) 280 PetscCallA(DMPlexGetPartitioner(dm,part,ierr)) 281 PetscCallA(PetscPartitionerSetFromOptions(part,ierr)) 282 PetscCallA(DMPlexDistribute(dm,0_kPI,migrationSF,pdm,ierr)) 283 284 if (numProc > 1) then 285 PetscCallA(DMPlexSetMigrationSF(pdm,migrationSF,ierr)) 286 PetscCallA(PetscSFDestroy(migrationSF,ierr)) 287 else 288 pdm = dm 289 end if 290 PetscCallA(DMViewFromOptions(pdm,PETSC_NULL_OPTIONS,'-dm_view',ierr)) 291 292 ! Get DM and IS for each field of dm 293 PetscCallA(DMCreateSubDM(pdm, 1_kPI, fieldU, isU, dmU,ierr)) 294 PetscCallA(DMCreateSubDM(pdm, 1_kPI, fieldA, isA, dmA,ierr)) 295 PetscCallA(DMCreateSubDM(pdm, 1_kPI, fieldS, isS, dmS,ierr)) 296 PetscCallA(DMCreateSubDM(pdm, 2_kPI, fieldUA, isUA, dmUA,ierr)) 297 298 !Create the exodus result file 299 allocate(dmList(2)) 300 dmList(1) = dmU; 301 dmList(2) = dmA; 302 PetscCallA(DMCreateSuperDM(dmList,2_kPI,PETSC_NULL_IS,dmUA2,ierr)) 303 deallocate(dmList) 304 305 PetscCallA(DMGetGlobalVector(pdm, X,ierr)) 306 PetscCallA(DMGetGlobalVector(dmU, U,ierr)) 307 PetscCallA(DMGetGlobalVector(dmA, A,ierr)) 308 PetscCallA(DMGetGlobalVector(dmS, S,ierr)) 309 PetscCallA(DMGetGlobalVector(dmUA, UA,ierr)) 310 PetscCallA(DMGetGlobalVector(dmUA2, UA2,ierr)) 311 312 PetscCallA(PetscObjectSetName(U, 'U',ierr)) 313 PetscCallA(PetscObjectSetName(A, 'Alpha',ierr)) 314 PetscCallA(PetscObjectSetName(S, 'Sigma',ierr)) 315 PetscCallA(PetscObjectSetName(UA, 'UAlpha',ierr)) 316 PetscCallA(PetscObjectSetName(UA2, 'UAlpha2',ierr)) 317 PetscCallA(VecSet(X, -111.0_kPR,ierr)) 318 319 ! 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 */ 320 PetscCallA(DMGetLocalSection(dmUA, sectionUA,ierr)) 321 PetscCallA(DMGetLocalVector(dmUA, UALoc,ierr)) 322 PetscCallA(VecGetArrayF90(UALoc, cval,ierr)) 323 PetscCallA(DMGetCoordinateSection(dmUA, coordSection,ierr)) 324 PetscCallA(DMGetCoordinatesLocal(dmUA, coord,ierr)) 325 PetscCallA(DMPlexGetChart(dmUA, pStart, pEnd,ierr)) 326 327 do p = pStart,pEnd-1 328 PetscCallA(PetscSectionGetDof(sectionUA, p, dofUA,ierr)) 329 if (dofUA > 0) then 330 PetscCallA(PetscSectionGetOffset(sectionUA, p, offUA,ierr)) 331 PetscCallA(DMPlexVecGetClosure(dmUA, coordSection, coord, p, xyz,ierr)) 332 closureSize = size(xyz) 333 do i = 1,sdim 334 do j = 0, closureSize-1,sdim 335 cval(offUA+i) = cval(offUA+i) + xyz(j/sdim+i) 336 end do 337 cval(offUA+i) = cval(offUA+i) * sdim / closureSize; 338 cval(offUA+sdim+1) = cval(offUA+sdim+1) + cval(offUA+i)**2 339 end do 340 PetscCallA(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, xyz,ierr)) 341 end if 342 end do 343 344 PetscCallA(VecRestoreArrayF90(UALoc, cval,ierr)) 345 PetscCallA(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA,ierr)) 346 PetscCallA(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA,ierr)) 347 PetscCallA(DMRestoreLocalVector(dmUA, UALoc,ierr)) 348 349 !Update X 350 PetscCallA(VecISCopy(X, isUA, SCATTER_FORWARD, UA,ierr)) 351 ! Restrict to U and Alpha 352 PetscCallA(VecISCopy(X, isU, SCATTER_REVERSE, U,ierr)) 353 PetscCallA(VecISCopy(X, isA, SCATTER_REVERSE, A,ierr)) 354 PetscCallA(VecViewFromOptions(UA, PETSC_NULL_OPTIONS, '-ua_vec_view',ierr)) 355 PetscCallA(VecViewFromOptions(U, PETSC_NULL_OPTIONS, '-u_vec_view',ierr)) 356 PetscCallA(VecViewFromOptions(A, PETSC_NULL_OPTIONS, '-a_vec_view',ierr)) 357 ! restrict to UA2 358 PetscCallA(VecISCopy(X, isUA, SCATTER_REVERSE, UA2,ierr)) 359 PetscCallA(VecViewFromOptions(UA2, PETSC_NULL_OPTIONS, '-ua2_vec_view',ierr)) 360 361 ! Getting Natural Vec 362 PetscCallA(DMSetOutputSequenceNumber(dmU, 0_kPI, time, ierr)) 363 PetscCallA(DMSetOutputSequenceNumber(dmA, 0_kPI, time, ierr)) 364 365 PetscCallA(VecView(U, viewer,ierr)) 366 PetscCallA(VecView(A, viewer,ierr)) 367 368 !Saving U and Alpha in one shot. 369 !For this, we need to cheat and change the Vec's name 370 !Note that in the end we write variables one component at a time, 371 !so that there is no real value in doing this 372 PetscCallA(DMSetOutputSequenceNumber(dmUA,1_kPI,time,ierr)) 373 PetscCallA(DMGetGlobalVector(dmUA, tmpVec,ierr)) 374 PetscCallA(VecCopy(UA, tmpVec,ierr)) 375 PetscCallA(PetscObjectSetName(tmpVec, 'U',ierr)) 376 PetscCallA(VecView(tmpVec, viewer,ierr)) 377 378 ! Reading nodal variables in Exodus file 379 PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr)) 380 PetscCallA(VecLoad(tmpVec, viewer,ierr)) 381 PetscCallA(VecAXPY(UA, -1.0_kPR, tmpVec,ierr)) 382 PetscCallA(VecNorm(UA, NORM_INFINITY, norm,ierr)) 383 if (norm > PETSC_SQRT_MACHINE_EPSILON) then 384 write(IOBuffer,'("UAlpha ||Vin - Vout|| = ",ES12.5)') norm 385 end if 386 PetscCallA(DMRestoreGlobalVector(dmUA, tmpVec,ierr)) 387 388 ! same thing with the UA2 Vec obtained from the superDM 389 PetscCallA(DMGetGlobalVector(dmUA2, tmpVec,ierr)) 390 PetscCallA(VecCopy(UA2, tmpVec,ierr)) 391 PetscCallA(PetscObjectSetName(tmpVec, 'U',ierr)) 392 PetscCallA(DMSetOutputSequenceNumber(dmUA2,2_kPI,time,ierr)) 393 PetscCallA(VecView(tmpVec, viewer,ierr)) 394 395 ! Reading nodal variables in Exodus file 396 PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr)) 397 PetscCallA(VecLoad(tmpVec,viewer,ierr)) 398 PetscCallA(VecAXPY(UA2, -1.0_kPR, tmpVec,ierr)) 399 PetscCallA(VecNorm(UA2, NORM_INFINITY, norm,ierr)) 400 if (norm > PETSC_SQRT_MACHINE_EPSILON) then 401 write(IOBuffer,'("UAlpha2 ||Vin - Vout|| = ",ES12.5)') norm 402 end if 403 PetscCallA(DMRestoreGlobalVector(dmUA2, tmpVec,ierr)) 404 405 ! Building and saving Sigma 406 ! We set sigma_0 = rank (to see partitioning) 407 ! sigma_1 = cell set ID 408 ! sigma_2 = x_coordinate of the cell center of mass 409 PetscCallA(DMGetCoordinateSection(dmS, coordSection,ierr)) 410 PetscCallA(DMGetCoordinatesLocal(dmS, coord,ierr)) 411 PetscCallA(DMGetLabelIdIS(dmS, 'Cell Sets', csIS,ierr)) 412 PetscCallA(DMGetLabelSize(dmS, 'Cell Sets',numCS,ierr)) 413 PetscCallA(ISGetIndicesF90(csIS, csID,ierr)) 414 415 do set = 1, numCS 416 PetscCallA(DMGetStratumIS(dmS, 'Cell Sets', csID(set), cellIS,ierr)) 417 PetscCallA(ISGetIndicesF90(cellIS, cellID,ierr)) 418 PetscCallA(ISGetSize(cellIS, numCells,ierr)) 419 do cell = 1,numCells 420 PetscCallA(DMPlexVecGetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval,ierr)) 421 PetscCallA(DMPlexVecGetClosure(dmS, coordSection, coord, cellID(cell), xyz,ierr)) 422 cval(1) = rank 423 cval(2) = csID(set) 424 cval(3) = 0.0_kPR 425 do c = 1, size(xyz),sdim 426 cval(3) = cval(3) + xyz(c) 427 end do 428 cval(3) = cval(3) * sdim / size(xyz) 429 PetscCallA(DMPlexVecSetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval, INSERT_ALL_VALUES,ierr)) 430 PetscCallA(DMPlexVecRestoreClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval,ierr)) 431 PetscCallA(DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID(cell), xyz,ierr)) 432 end do 433 PetscCallA(ISRestoreIndicesF90(cellIS, cellID,ierr)) 434 PetscCallA(ISDestroy(cellIS,ierr)) 435 end do 436 PetscCallA(ISRestoreIndicesF90(csIS, csID,ierr)) 437 PetscCallA(ISDestroy(csIS,ierr)) 438 PetscCallA(VecViewFromOptions(S, PETSC_NULL_OPTIONS, '-s_vec_view',ierr)) 439 440 ! Writing zonal variables in Exodus file 441 PetscCallA(DMSetOutputSequenceNumber(dmS,0_kPI,time,ierr)) 442 PetscCallA(VecView(S,viewer,ierr)) 443 444 ! Reading zonal variables in Exodus file */ 445 PetscCallA(DMGetGlobalVector(dmS, tmpVec,ierr)) 446 PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr)) 447 PetscCallA(PetscObjectSetName(tmpVec, 'Sigma',ierr)) 448 PetscCallA(VecLoad(tmpVec,viewer,ierr)) 449 PetscCallA(VecAXPY(S, -1.0_kPR, tmpVec,ierr)) 450 PetscCallA(VecNorm(S, NORM_INFINITY,norm,ierr)) 451 if (norm > PETSC_SQRT_MACHINE_EPSILON) then 452 write(IOBuffer,'("Sigma ||Vin - Vout|| = ",ES12.5)') norm 453 end if 454 PetscCallA(DMRestoreGlobalVector(dmS, tmpVec,ierr)) 455 456 PetscCallA(DMRestoreGlobalVector(dmUA2, UA2,ierr)) 457 PetscCallA(DMRestoreGlobalVector(dmUA, UA,ierr)) 458 PetscCallA(DMRestoreGlobalVector(dmS, S,ierr)) 459 PetscCallA(DMRestoreGlobalVector(dmA, A,ierr)) 460 PetscCallA(DMRestoreGlobalVector(dmU, U,ierr)) 461 PetscCallA(DMRestoreGlobalVector(pdm, X,ierr)) 462 PetscCallA(DMDestroy(dmU,ierr)) 463 PetscCallA(ISDestroy(isU,ierr)) 464 PetscCallA(DMDestroy(dmA,ierr)) 465 PetscCallA(ISDestroy(isA,ierr)) 466 PetscCallA(DMDestroy(dmS,ierr)) 467 PetscCallA(ISDestroy(isS,ierr)) 468 PetscCallA(DMDestroy(dmUA,ierr)) 469 PetscCallA(ISDestroy(isUA,ierr)) 470 PetscCallA(DMDestroy(dmUA2,ierr)) 471 PetscCallA(DMDestroy(pdm,ierr)) 472 if (numProc > 1) then 473 PetscCallA(DMDestroy(dm,ierr)) 474 end if 475 476 deallocate(pStartDepth) 477 deallocate(pEndDepth) 478 479 PetscCallA(PetscViewerDestroy(viewer,ierr)) 480 PetscCallA(PetscFinalize(ierr)) 481end program ex62f90 482 483! /*TEST 484! 485! build: 486! requires: exodusii pnetcdf !complex 487! # 2D seq 488! test: 489! suffix: 0 490! 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 491! #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 492! test: 493! suffix: 1 494! 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 495! 496! test: 497! suffix: 2 498! 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 499! #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 500! test: 501! suffix: 3 502! 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 503! test: 504! suffix: 4 505! 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 506! test: 507! suffix: 5 508! 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 509! # 2D par 510! test: 511! suffix: 6 512! nsize: 2 513! 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 514! #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 515! test: 516! suffix: 7 517! nsize: 2 518! 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 519! test: 520! suffix: 8 521! nsize: 2 522! 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 523! #TODO: bug in PetscCallA(to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name 524! test: 525! suffix: 9 526! nsize: 2 527! 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 528! test: 529! suffix: 10 530! nsize: 2 531! 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 532! test: 533! # Something is now broken with parallel read/write for whatever shape H is 534! TODO: broken 535! suffix: 11 536! nsize: 2 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 539! #3d seq 540! test: 541! suffix: 12 542! 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 543! test: 544! suffix: 13 545! 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 546! #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 547! test: 548! suffix: 14 549! 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 550! test: 551! suffix: 15 552! 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 553! #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 554! #3d par 555! test: 556! suffix: 16 557! nsize: 2 558! 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 559! test: 560! suffix: 17 561! nsize: 2 562! 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 563! #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 564! test: 565! suffix: 18 566! nsize: 2 567! 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 568! test: 569! suffix: 19 570! nsize: 2 571! 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 572! #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 573 574! TEST*/ 575