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