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 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 PetscObjectGetComm(dm,comm,ierr);CHKERRA(ierr) 171 call PetscSectionCreate(comm, section,ierr);CHKERRA(ierr) 172 call PetscSectionSetNumFields(section, 3_kPI,ierr);CHKERRA(ierr) 173 call PetscSectionSetFieldName(section, fieldU, "U",ierr);CHKERRA(ierr) 174 call PetscSectionSetFieldName(section, fieldA, "Alpha",ierr);CHKERRA(ierr) 175 call PetscSectionSetFieldName(section, fieldS, "Sigma",ierr);CHKERRA(ierr) 176 call DMPlexGetChart(dm, pStart, pEnd,ierr);CHKERRA(ierr) 177 call PetscSectionSetChart(section, pStart, pEnd,ierr);CHKERRA(ierr) 178 179 allocate(pStartDepth(sdim+1)) 180 allocate(pEndDepth(sdim+1)) 181 do d = 1, sdim+1 182 call DMPlexGetDepthStratum(dm, d-1, pStartDepth(d), pEndDepth(d),ierr);CHKERRA(ierr) 183 end do 184 185 ! Vector field U, Scalar field Alpha, Tensor field Sigma 186 call PetscSectionSetFieldComponents(section, fieldU, sdim,ierr);CHKERRA(ierr); 187 call PetscSectionSetFieldComponents(section, fieldA, 1_kPI,ierr);CHKERRA(ierr); 188 call PetscSectionSetFieldComponents(section, fieldS, sdim*(sdim+1)/2,ierr);CHKERRA(ierr); 189 190 ! Going through cell sets then cells, and setting up storage for the sections 191 call DMGetLabelSize(dm, "Cell Sets", numCS, ierr);CHKERRA(ierr) 192 call DMGetLabelIdIS(dm, "Cell Sets", csIS, ierr);CHKERRA(ierr) 193 call ISGetIndicesF90(csIS, csID, ierr);CHKERRA(ierr) 194 do set = 1,numCS 195 call DMGetStratumSize(dm, "Cell Sets", csID(set), numCells,ierr);CHKERRA(ierr) 196 call DMGetStratumIS(dm, "Cell Sets", csID(set), cellIS,ierr);CHKERRA(ierr) 197 if (numCells > 0) then 198 select case(sdim) 199 case(2) 200 dofs => dofS2D 201 case(3) 202 dofs => dofS3D 203 case default 204 write(IOBuffer,'("No layout for dimension ",I2)') sdim 205 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,IOBuffer) 206 end select ! sdim 207 208 ! Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes 209 ! It will not be enough to identify more exotic elements like pyramid or prisms... */ 210 call ISGetIndicesF90(cellIS, cellID,ierr);CHKERRA(ierr) 211 nullify(closureA) 212 call DMPlexGetTransitiveClosure(dm,cellID(1), PETSC_TRUE, closureA,ierr);CHKERRA(ierr) 213 select case(size(closureA)/2) 214 case(7) ! Tri 215 if (order == 1) then 216 dofU => dofUP1Tri 217 dofA => dofAP1Tri 218 else 219 dofU => dofUP2Tri 220 dofA => dofAP2Tri 221 end if 222 case(9) ! Quad 223 if (order == 1) then 224 dofU => dofUP1Quad 225 dofA => dofAP1Quad 226 else 227 dofU => dofUP2Quad 228 dofA => dofAP2Quad 229 end if 230 case(15) ! Tet 231 if (order == 1) then 232 dofU => dofUP1Tet 233 dofA => dofAP1Tet 234 else 235 dofU => dofUP2Tet 236 dofA => dofAP2Tet 237 end if 238 case(27) ! Hex 239 if (order == 1) then 240 dofU => dofUP1Hex 241 dofA => dofAP1Hex 242 else 243 dofU => dofUP2Hex 244 dofA => dofAP2Hex 245 end if 246 case default 247 write(IOBuffer,'("Unknown element with closure size ",I2)') size(closureA)/2 248 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,IOBuffer) 249 end select 250 call DMPlexRestoreTransitiveClosure(dm, cellID(1), PETSC_TRUE,closureA,ierr);CHKERRA(ierr) 251 do cell = 1,numCells! 252 nullify(closure) 253 call DMPlexGetTransitiveClosure(dm, cellID(cell), PETSC_TRUE, closure,ierr);CHKERRA(ierr) 254 do p = 1,size(closure),2 255 ! find the depth of p 256 do d = 1,sdim+1 257 if ((closure(p) >= pStartDepth(d)) .and. (closure(p) < pEndDepth(d))) then 258 call PetscSectionSetDof(section, closure(p), dofU(d)+dofA(d)+dofS(d),ierr);CHKERRA(ierr) 259 call PetscSectionSetFieldDof(section, closure(p), fieldU, dofU(d),ierr);CHKERRA(ierr) 260 call PetscSectionSetFieldDof(section, closure(p), fieldA, dofA(d),ierr);CHKERRA(ierr) 261 call PetscSectionSetFieldDof(section, closure(p), fieldS, dofS(d),ierr);CHKERRA(ierr) 262 end if ! closure(p) 263 end do ! d 264 end do ! p 265 call DMPlexRestoreTransitiveClosure(dm, cellID(cell), PETSC_TRUE, closure,ierr);CHKERRA(ierr) 266 end do ! cell 267 call ISRestoreIndicesF90(cellIS, cellID,ierr);CHKERRA(ierr) 268 call ISDestroy(cellIS,ierr);CHKERRA(ierr) 269 end if ! numCells 270 end do ! set 271 call ISRestoreIndicesF90(csIS, csID,ierr);CHKERRA(ierr) 272 call ISDestroy(csIS,ierr);CHKERRA(ierr) 273 call PetscSectionSetUp(section,ierr);CHKERRA(ierr) 274 call DMSetLocalSection(dm, section,ierr);CHKERRA(ierr) 275 call PetscObjectViewFromOptions(section, PETSC_NULL_SECTION, "-dm_section_view",ierr);CHKERRQ(ierr) 276 call PetscSectionDestroy(section,ierr);CHKERRA(ierr) 277 278 call DMSetUseNatural(dm,PETSC_TRUE,ierr);CHKERRA(ierr) 279 call DMPlexGetPartitioner(dm,part,ierr);CHKERRA(ierr) 280 call PetscPartitionerSetFromOptions(part,ierr);CHKERRA(ierr) 281 call DMPlexDistribute(dm,0_kPI,migrationSF,pdm,ierr);CHKERRA(ierr) 282 283 if (numProc > 1) then 284 call DMPlexSetMigrationSF(pdm,migrationSF,ierr);CHKERRA(ierr) 285 call PetscSFDestroy(migrationSF,ierr);CHKERRA(ierr) 286 call DMDestroy(dm,ierr);CHKERRA(ierr) 287 dm = pdm 288 end if 289 call DMViewFromOptions(dm,PETSC_NULL_OPTIONS,"-dm_view",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 call DMSetOutputSequenceNumber(dmUA,1_kPI,time,ierr);CHKERRA(ierr) 372 call DMGetGlobalVector(dmUA, tmpVec,ierr);CHKERRA(ierr) 373 call VecCopy(UA, tmpVec,ierr);CHKERRA(ierr) 374 call PetscObjectSetName(tmpVec, "U",ierr);CHKERRA(ierr) 375 call VecView(tmpVec, viewer,ierr);CHKERRA(ierr) 376 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 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 ex26f90 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