1program ex26f90 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(ierr)) 81 82 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)) 83 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,numProc,ierr)) 84 PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-i",ifilename,flg,ierr)) 85 if (.not. flg) then 86 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"missing input file name -i <input file name>") 87 end if 88 PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-o",ofilename,flg,ierr)) 89 if (.not. flg) then 90 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"missing output file name -o <output file name>") 91 end if 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 SETERRQ(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 PetscCallA(DMGetStratumSize(dm, "Cell Sets", csID(set), numCells,ierr)) 193 PetscCallA(DMGetStratumIS(dm, "Cell Sets", csID(set), cellIS,ierr)) 194 if (numCells > 0) then 195 select case(sdim) 196 case(2) 197 dofs => dofS2D 198 case(3) 199 dofs => dofS3D 200 case default 201 write(IOBuffer,'("No layout for dimension ",I2)') sdim 202 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,IOBuffer) 203 end select ! sdim 204 205 ! Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes 206 ! It will not be enough to identify more exotic elements like pyramid or prisms... */ 207 PetscCallA(ISGetIndicesF90(cellIS, cellID,ierr)) 208 nullify(closureA) 209 PetscCallA(DMPlexGetTransitiveClosure(dm,cellID(1), PETSC_TRUE, closureA,ierr)) 210 select case(size(closureA)/2) 211 case(7) ! Tri 212 if (order == 1) then 213 dofU => dofUP1Tri 214 dofA => dofAP1Tri 215 else 216 dofU => dofUP2Tri 217 dofA => dofAP2Tri 218 end if 219 case(9) ! Quad 220 if (order == 1) then 221 dofU => dofUP1Quad 222 dofA => dofAP1Quad 223 else 224 dofU => dofUP2Quad 225 dofA => dofAP2Quad 226 end if 227 case(15) ! Tet 228 if (order == 1) then 229 dofU => dofUP1Tet 230 dofA => dofAP1Tet 231 else 232 dofU => dofUP2Tet 233 dofA => dofAP2Tet 234 end if 235 case(27) ! Hex 236 if (order == 1) then 237 dofU => dofUP1Hex 238 dofA => dofAP1Hex 239 else 240 dofU => dofUP2Hex 241 dofA => dofAP2Hex 242 end if 243 case default 244 write(IOBuffer,'("Unknown element with closure size ",I2)') size(closureA)/2 245 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,IOBuffer) 246 end select 247 PetscCallA(DMPlexRestoreTransitiveClosure(dm, cellID(1), PETSC_TRUE,closureA,ierr)) 248 do cell = 1,numCells! 249 nullify(closure) 250 PetscCallA(DMPlexGetTransitiveClosure(dm, cellID(cell), PETSC_TRUE, closure,ierr)) 251 do p = 1,size(closure),2 252 ! find the depth of p 253 do d = 1,sdim+1 254 if ((closure(p) >= pStartDepth(d)) .and. (closure(p) < pEndDepth(d))) then 255 PetscCallA(PetscSectionSetDof(section, closure(p), dofU(d)+dofA(d)+dofS(d),ierr)) 256 PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldU, dofU(d),ierr)) 257 PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldA, dofA(d),ierr)) 258 PetscCallA(PetscSectionSetFieldDof(section, closure(p), fieldS, dofS(d),ierr)) 259 end if ! closure(p) 260 end do ! d 261 end do ! p 262 PetscCallA(DMPlexRestoreTransitiveClosure(dm, cellID(cell), PETSC_TRUE, closure,ierr)) 263 end do ! cell 264 PetscCallA(ISRestoreIndicesF90(cellIS, cellID,ierr)) 265 PetscCallA(ISDestroy(cellIS,ierr)) 266 end if ! numCells 267 end do ! set 268 PetscCallA(ISRestoreIndicesF90(csIS, csID,ierr)) 269 PetscCallA(ISDestroy(csIS,ierr)) 270 PetscCallA(PetscSectionSetUp(section,ierr)) 271 PetscCallA(DMSetLocalSection(dm, section,ierr)) 272 PetscCallA(PetscObjectViewFromOptions(section, PETSC_NULL_SECTION, "-dm_section_view",ierr)) 273 PetscCallA(PetscSectionDestroy(section,ierr)) 274 275 PetscCallA(DMSetUseNatural(dm,PETSC_TRUE,ierr)) 276 PetscCallA(DMPlexGetPartitioner(dm,part,ierr)) 277 PetscCallA(PetscPartitionerSetFromOptions(part,ierr)) 278 PetscCallA(DMPlexDistribute(dm,0_kPI,migrationSF,pdm,ierr)) 279 280 if (numProc > 1) then 281 PetscCallA(DMPlexSetMigrationSF(pdm,migrationSF,ierr)) 282 PetscCallA(PetscSFDestroy(migrationSF,ierr)) 283 PetscCallA(DMDestroy(dm,ierr)) 284 dm = pdm 285 end if 286 PetscCallA(DMViewFromOptions(dm,PETSC_NULL_OPTIONS,"-dm_view",ierr)) 287 288 ! Get DM and IS for each field of dm 289 PetscCallA(DMCreateSubDM(dm, 1_kPI, fieldU, isU, dmU,ierr)) 290 PetscCallA(DMCreateSubDM(dm, 1_kPI, fieldA, isA, dmA,ierr)) 291 PetscCallA(DMCreateSubDM(dm, 1_kPI, fieldS, isS, dmS,ierr)) 292 PetscCallA(DMCreateSubDM(dm, 2_kPI, fieldUA, isUA, dmUA,ierr)) 293 294 !Create the exodus result file 295 allocate(dmList(2)) 296 dmList(1) = dmU; 297 dmList(2) = dmA; 298 PetscCallA(DMCreateSuperDM(dmList,2_kPI,PETSC_NULL_IS,dmUA2,ierr)) 299 deallocate(dmList) 300 301 PetscCallA(DMGetGlobalVector(dm, X,ierr)) 302 PetscCallA(DMGetGlobalVector(dmU, U,ierr)) 303 PetscCallA(DMGetGlobalVector(dmA, A,ierr)) 304 PetscCallA(DMGetGlobalVector(dmS, S,ierr)) 305 PetscCallA(DMGetGlobalVector(dmUA, UA,ierr)) 306 PetscCallA(DMGetGlobalVector(dmUA2, UA2,ierr)) 307 308 PetscCallA(PetscObjectSetName(U, "U",ierr)) 309 PetscCallA(PetscObjectSetName(A, "Alpha",ierr)) 310 PetscCallA(PetscObjectSetName(S, "Sigma",ierr)) 311 PetscCallA(PetscObjectSetName(UA, "UAlpha",ierr)) 312 PetscCallA(PetscObjectSetName(UA2, "UAlpha2",ierr)) 313 PetscCallA(VecSet(X, -111.0_kPR,ierr)) 314 315 ! 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 */ 316 PetscCallA(DMGetLocalSection(dmUA, sectionUA,ierr)) 317 PetscCallA(DMGetLocalVector(dmUA, UALoc,ierr)) 318 PetscCallA(VecGetArrayF90(UALoc, cval,ierr)) 319 PetscCallA(DMGetCoordinateSection(dmUA, coordSection,ierr)) 320 PetscCallA(DMGetCoordinatesLocal(dmUA, coord,ierr)) 321 PetscCallA(DMPlexGetChart(dmUA, pStart, pEnd,ierr)) 322 323 do p = pStart,pEnd-1 324 PetscCallA(PetscSectionGetDof(sectionUA, p, dofUA,ierr)) 325 if (dofUA > 0) then 326 PetscCallA(PetscSectionGetOffset(sectionUA, p, offUA,ierr)) 327 PetscCallA(DMPlexVecGetClosure(dmUA, coordSection, coord, p, xyz,ierr)) 328 closureSize = size(xyz) 329 do i = 1,sdim 330 do j = 0, closureSize-1,sdim 331 cval(offUA+i) = cval(offUA+i) + xyz(j/sdim+i) 332 end do 333 cval(offUA+i) = cval(offUA+i) * sdim / closureSize; 334 cval(offUA+sdim+1) = cval(offUA+sdim+1) + cval(offUA+i)**2 335 end do 336 PetscCallA(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, xyz,ierr)) 337 end if 338 end do 339 340 PetscCallA(VecRestoreArrayF90(UALoc, cval,ierr)) 341 PetscCallA(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA,ierr)) 342 PetscCallA(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA,ierr)) 343 PetscCallA(DMRestoreLocalVector(dmUA, UALoc,ierr)) 344 345 !Update X 346 PetscCallA(VecISCopy(X, isUA, SCATTER_FORWARD, UA,ierr)) 347 ! Restrict to U and Alpha 348 PetscCallA(VecISCopy(X, isU, SCATTER_REVERSE, U,ierr)) 349 PetscCallA(VecISCopy(X, isA, SCATTER_REVERSE, A,ierr)) 350 PetscCallA(VecViewFromOptions(UA, PETSC_NULL_OPTIONS, "-ua_vec_view",ierr)) 351 PetscCallA(VecViewFromOptions(U, PETSC_NULL_OPTIONS, "-u_vec_view",ierr)) 352 PetscCallA(VecViewFromOptions(A, PETSC_NULL_OPTIONS, "-a_vec_view",ierr)) 353 ! restrict to UA2 354 PetscCallA(VecISCopy(X, isUA, SCATTER_REVERSE, UA2,ierr)) 355 PetscCallA(VecViewFromOptions(UA2, PETSC_NULL_OPTIONS, "-ua2_vec_view",ierr)) 356 357 ! Writing nodal variables to ExodusII file 358 PetscCallA(DMSetOutputSequenceNumber(dmU,0_kPI,time,ierr)) 359 PetscCallA(DMSetOutputSequenceNumber(dmA,0_kPI,time,ierr)) 360 361 PetscCallA(VecView(U, viewer,ierr)) 362 PetscCallA(VecView(A, viewer,ierr)) 363 364 ! Saving U and Alpha in one shot. 365 ! For this, we need to cheat and change the Vec's name 366 ! Note that in the end we write variables one component at a time, 367 ! so that there is no real value in doing this 368 PetscCallA(DMSetOutputSequenceNumber(dmUA,1_kPI,time,ierr)) 369 PetscCallA(DMGetGlobalVector(dmUA, tmpVec,ierr)) 370 PetscCallA(VecCopy(UA, tmpVec,ierr)) 371 PetscCallA(PetscObjectSetName(tmpVec, "U",ierr)) 372 PetscCallA(VecView(tmpVec, viewer,ierr)) 373 374 ! Reading nodal variables in Exodus file 375 PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr)) 376 PetscCallA(VecLoad(tmpVec, viewer,ierr)) 377 PetscCallA(VecAXPY(UA, -1.0_kPR, tmpVec,ierr)) 378 PetscCallA(VecNorm(UA, NORM_INFINITY, norm,ierr)) 379 if (norm > PETSC_SQRT_MACHINE_EPSILON) then 380 write(IOBuffer,'("UAlpha ||Vin - Vout|| = ",ES12.5)') norm 381 end if 382 PetscCallA(DMRestoreGlobalVector(dmUA, tmpVec,ierr)) 383 384 ! same thing with the UA2 Vec obtained from the superDM 385 PetscCallA(DMGetGlobalVector(dmUA2, tmpVec,ierr)) 386 PetscCallA(VecCopy(UA2, tmpVec,ierr)) 387 PetscCallA(PetscObjectSetName(tmpVec, "U",ierr)) 388 PetscCallA(DMSetOutputSequenceNumber(dmUA2,2_kPI,time,ierr)) 389 PetscCallA(VecView(tmpVec, viewer,ierr)) 390 391 ! Reading nodal variables in Exodus file 392 PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr)) 393 PetscCallA(VecLoad(tmpVec,viewer,ierr)) 394 PetscCallA(VecAXPY(UA2, -1.0_kPR, tmpVec,ierr)) 395 PetscCallA(VecNorm(UA2, NORM_INFINITY, norm,ierr)) 396 if (norm > PETSC_SQRT_MACHINE_EPSILON) then 397 write(IOBuffer,'("UAlpha2 ||Vin - Vout|| = ",ES12.5)') norm 398 end if 399 PetscCallA(DMRestoreGlobalVector(dmUA2, tmpVec,ierr)) 400 401 ! Building and saving Sigma 402 ! We set sigma_0 = rank (to see partitioning) 403 ! sigma_1 = cell set ID 404 ! sigma_2 = x_coordinate of the cell center of mass 405 PetscCallA(DMGetCoordinateSection(dmS, coordSection,ierr)) 406 PetscCallA(DMGetCoordinatesLocal(dmS, coord,ierr)) 407 PetscCallA(DMGetLabelIdIS(dmS, "Cell Sets", csIS,ierr)) 408 PetscCallA(DMGetLabelSize(dmS, "Cell Sets",numCS,ierr)) 409 PetscCallA(ISGetIndicesF90(csIS, csID,ierr)) 410 411 do set = 1, numCS 412 PetscCallA(DMGetStratumIS(dmS, "Cell Sets", csID(set), cellIS,ierr)) 413 PetscCallA(ISGetIndicesF90(cellIS, cellID,ierr)) 414 PetscCallA(ISGetSize(cellIS, numCells,ierr)) 415 do cell = 1,numCells 416 PetscCallA(DMPlexVecGetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval,ierr)) 417 PetscCallA(DMPlexVecGetClosure(dmS, coordSection, coord, cellID(cell), xyz,ierr)) 418 cval(1) = rank 419 cval(2) = csID(set) 420 cval(3) = 0.0_kPR 421 do c = 1, size(xyz),sdim 422 cval(3) = cval(3) + xyz(c) 423 end do 424 cval(3) = cval(3) * sdim / size(xyz) 425 PetscCallA(DMPlexVecSetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval, INSERT_ALL_VALUES,ierr)) 426 PetscCallA(DMPlexVecRestoreClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval,ierr)) 427 PetscCallA(DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID(cell), xyz,ierr)) 428 end do 429 PetscCallA(ISRestoreIndicesF90(cellIS, cellID,ierr)) 430 PetscCallA(ISDestroy(cellIS,ierr)) 431 end do 432 PetscCallA(ISRestoreIndicesF90(csIS, csID,ierr)) 433 PetscCallA(ISDestroy(csIS,ierr)) 434 PetscCallA(VecViewFromOptions(S, PETSC_NULL_OPTIONS, "-s_vec_view",ierr)) 435 436 ! Writing zonal variables in Exodus file 437 PetscCallA(DMSetOutputSequenceNumber(dmS,0_kPI,time,ierr)) 438 PetscCallA(VecView(S,viewer,ierr)) 439 440 ! Reading zonal variables in Exodus file */ 441 PetscCallA(DMGetGlobalVector(dmS, tmpVec,ierr)) 442 PetscCallA(VecSet(tmpVec, -1000.0_kPR,ierr)) 443 PetscCallA(PetscObjectSetName(tmpVec, "Sigma",ierr)) 444 PetscCallA(VecLoad(tmpVec,viewer,ierr)) 445 PetscCallA(VecAXPY(S, -1.0_kPR, tmpVec,ierr)) 446 PetscCallA(VecNorm(S, NORM_INFINITY,norm,ierr)) 447 if (norm > PETSC_SQRT_MACHINE_EPSILON) then 448 write(IOBuffer,'("Sigma ||Vin - Vout|| = ",ES12.5)') norm 449 end if 450 PetscCallA(DMRestoreGlobalVector(dmS, tmpVec,ierr)) 451 452 PetscCallA(DMRestoreGlobalVector(dmUA2, UA2,ierr)) 453 PetscCallA(DMRestoreGlobalVector(dmUA, UA,ierr)) 454 PetscCallA(DMRestoreGlobalVector(dmS, S,ierr)) 455 PetscCallA(DMRestoreGlobalVector(dmA, A,ierr)) 456 PetscCallA(DMRestoreGlobalVector(dmU, U,ierr)) 457 PetscCallA(DMRestoreGlobalVector(dm, X,ierr)) 458 PetscCallA(DMDestroy(dmU,ierr)) 459 PetscCallA(ISDestroy(isU,ierr)) 460 PetscCallA(DMDestroy(dmA,ierr)) 461 PetscCallA(ISDestroy(isA,ierr)) 462 PetscCallA(DMDestroy(dmS,ierr)) 463 PetscCallA(ISDestroy(isS,ierr)) 464 PetscCallA(DMDestroy(dmUA,ierr)) 465 PetscCallA(ISDestroy(isUA,ierr)) 466 PetscCallA(DMDestroy(dmUA2,ierr)) 467 PetscCallA(DMDestroy(dm,ierr)) 468 469 deallocate(pStartDepth) 470 deallocate(pEndDepth) 471 472 PetscCallA(PetscViewerDestroy(viewer,ierr)) 473 PetscCallA(PetscFinalize(ierr)) 474end program ex26f90 475 476! /*TEST 477! 478! build: 479! requires: exodusii pnetcdf !complex 480! # 2D seq 481! test: 482! suffix: 0 483! 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 484! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 485! test: 486! suffix: 1 487! 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 488! 489! test: 490! suffix: 2 491! 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 492! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 493! test: 494! suffix: 3 495! 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 496! test: 497! suffix: 4 498! 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 499! test: 500! suffix: 5 501! 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 502! # 2D par 503! test: 504! suffix: 6 505! nsize: 2 506! 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 507! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 508! test: 509! suffix: 7 510! nsize: 2 511! 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 512! test: 513! suffix: 8 514! nsize: 2 515! 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 516! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name 517! test: 518! suffix: 9 519! nsize: 2 520! 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 521! test: 522! suffix: 10 523! nsize: 2 524! 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 525! test: 526! # Something is now broken with parallel read/write for wahtever shape H is 527! TODO: broken 528! suffix: 11 529! nsize: 2 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 532! #3d seq 533! test: 534! suffix: 12 535! 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 536! test: 537! suffix: 13 538! 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 539! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 540! test: 541! suffix: 14 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 2 543! test: 544! suffix: 15 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 2 546! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 547! #3d par 548! test: 549! suffix: 16 550! nsize: 2 551! 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 552! test: 553! suffix: 17 554! nsize: 2 555! 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 556! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 557! test: 558! suffix: 18 559! nsize: 2 560! 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 561! test: 562! suffix: 19 563! nsize: 2 564! 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 565! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: One or more variable sizes violate format constraints 566! TEST*/ 567