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