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_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 PetscObjectGetComm(dm,comm,ierr);CHKERRA(ierr) 172 call PetscSectionCreate(comm, section,ierr);CHKERRA(ierr) 173 call PetscSectionSetNumFields(section, 3_kPI,ierr);CHKERRA(ierr) 174 call PetscSectionSetFieldName(section, fieldU, "U",ierr);CHKERRA(ierr) 175 call PetscSectionSetFieldName(section, fieldA, "Alpha",ierr);CHKERRA(ierr) 176 call PetscSectionSetFieldName(section, fieldS, "Sigma",ierr);CHKERRA(ierr) 177 call DMPlexGetChart(dm, pStart, pEnd,ierr);CHKERRA(ierr) 178 call PetscSectionSetChart(section, pStart, pEnd,ierr);CHKERRA(ierr) 179 180 allocate(pStartDepth(sdim+1)) 181 allocate(pEndDepth(sdim+1)) 182 do d = 1, sdim+1 183 call DMPlexGetDepthStratum(dm, d-1, pStartDepth(d), pEndDepth(d),ierr);CHKERRA(ierr) 184 end do 185 186 ! Vector field U, Scalar field Alpha, Tensor field Sigma 187 call PetscSectionSetFieldComponents(section, fieldU, sdim,ierr);CHKERRA(ierr); 188 call PetscSectionSetFieldComponents(section, fieldA, 1_kPI,ierr);CHKERRA(ierr); 189 call PetscSectionSetFieldComponents(section, fieldS, sdim*(sdim+1)/2,ierr);CHKERRA(ierr); 190 191 ! Going through cell sets then cells, and setting up storage for the sections 192 call DMGetLabelSize(dm, "Cell Sets", numCS, ierr);CHKERRA(ierr) 193 call DMGetLabelIdIS(dm, "Cell Sets", csIS, ierr);CHKERRA(ierr) 194 call ISGetIndicesF90(csIS, csID, ierr);CHKERRA(ierr) 195 do set = 1,numCS 196 call DMGetStratumSize(dm, "Cell Sets", csID(set), numCells,ierr);CHKERRA(ierr) 197 call DMGetStratumIS(dm, "Cell Sets", csID(set), cellIS,ierr);CHKERRA(ierr) 198 if (numCells > 0) then 199 select case(sdim) 200 case(2) 201 dofs => dofS2D 202 case(3) 203 dofs => dofS3D 204 case default 205 write(IOBuffer,'("No layout for dimension ",I2)') sdim 206 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,IOBuffer) 207 end select ! sdim 208 209 ! Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes 210 ! It will not be enough to identify more exotic elements like pyramid or prisms... */ 211 call ISGetIndicesF90(cellIS, cellID,ierr);CHKERRA(ierr) 212 nullify(closureA) 213 call DMPlexGetTransitiveClosure(dm,cellID(1), PETSC_TRUE, closureA,ierr);CHKERRA(ierr) 214 select case(size(closureA)/2) 215 case(7) ! Tri 216 if (order == 1) then 217 dofU => dofUP1Tri 218 dofA => dofAP1Tri 219 else 220 dofU => dofUP2Tri 221 dofA => dofAP2Tri 222 end if 223 case(9) ! Quad 224 if (order == 1) then 225 dofU => dofUP1Quad 226 dofA => dofAP1Quad 227 else 228 dofU => dofUP2Quad 229 dofA => dofAP2Quad 230 end if 231 case(15) ! Tet 232 if (order == 1) then 233 dofU => dofUP1Tet 234 dofA => dofAP1Tet 235 else 236 dofU => dofUP2Tet 237 dofA => dofAP2Tet 238 end if 239 case(27) ! Hex 240 if (order == 1) then 241 dofU => dofUP1Hex 242 dofA => dofAP1Hex 243 else 244 dofU => dofUP2Hex 245 dofA => dofAP2Hex 246 end if 247 case default 248 write(IOBuffer,'("Unknown element with closure size ",I2)') size(closureA)/2 249 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,IOBuffer) 250 end select 251 call DMPlexRestoreTransitiveClosure(dm, cellID(1), PETSC_TRUE,closureA,ierr);CHKERRA(ierr) 252 do cell = 1,numCells! 253 nullify(closure) 254 call DMPlexGetTransitiveClosure(dm, cellID(cell), PETSC_TRUE, closure,ierr);CHKERRA(ierr) 255 do p = 1,size(closure),2 256 ! find the depth of p 257 do d = 1,sdim+1 258 if ((closure(p) >= pStartDepth(d)) .and. (closure(p) < pEndDepth(d))) then 259 call PetscSectionSetDof(section, closure(p), dofU(d)+dofA(d)+dofS(d),ierr);CHKERRA(ierr) 260 call PetscSectionSetFieldDof(section, closure(p), fieldU, dofU(d),ierr);CHKERRA(ierr) 261 call PetscSectionSetFieldDof(section, closure(p), fieldA, dofA(d),ierr);CHKERRA(ierr) 262 call PetscSectionSetFieldDof(section, closure(p), fieldS, dofS(d),ierr);CHKERRA(ierr) 263 end if ! closure(p) 264 end do ! d 265 end do ! p 266 call DMPlexRestoreTransitiveClosure(dm, cellID(cell), PETSC_TRUE, closure,ierr);CHKERRA(ierr) 267 end do ! cell 268 call ISRestoreIndicesF90(cellIS, cellID,ierr);CHKERRA(ierr) 269 call ISDestroy(cellIS,ierr);CHKERRA(ierr) 270 end if ! numCells 271 end do ! set 272 call ISRestoreIndicesF90(csIS, csID,ierr);CHKERRA(ierr) 273 call ISDestroy(csIS,ierr);CHKERRA(ierr) 274 call PetscSectionSetUp(section,ierr);CHKERRA(ierr) 275 call DMSetLocalSection(dm, section,ierr);CHKERRA(ierr) 276 call PetscObjectViewFromOptions(section, PETSC_NULL_SECTION, "-dm_section_view",ierr);CHKERRQ(ierr) 277 call PetscSectionDestroy(section,ierr);CHKERRA(ierr) 278 279 call DMSetUseNatural(dm,PETSC_TRUE,ierr);CHKERRA(ierr) 280 call DMPlexGetPartitioner(dm,part,ierr);CHKERRA(ierr) 281 call PetscPartitionerSetFromOptions(part,ierr);CHKERRA(ierr) 282 call DMPlexDistribute(dm,0_kPI,migrationSF,pdm,ierr);CHKERRA(ierr) 283 284 if (numProc > 1) then 285 call DMPlexSetMigrationSF(pdm,migrationSF,ierr);CHKERRA(ierr) 286 call PetscSFDestroy(migrationSF,ierr);CHKERRA(ierr) 287 call DMDestroy(dm,ierr);CHKERRA(ierr) 288 dm = pdm 289 end if 290 call DMViewFromOptions(dm,PETSC_NULL_OPTIONS,"-dm_view",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 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 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 395 ! Reading nodal variables in Exodus file 396 call VecSet(tmpVec, -1000.0_kPR,ierr);CHKERRA(ierr) 397 call VecLoad(tmpVec,viewer,ierr);CHKERRA(ierr) 398 call VecAXPY(UA2, -1.0_kPR, tmpVec,ierr);CHKERRA(ierr) 399 call VecNorm(UA2, NORM_INFINITY, norm,ierr);CHKERRA(ierr) 400 if (norm > PETSC_SQRT_MACHINE_EPSILON) then 401 write(IOBuffer,'("UAlpha2 ||Vin - Vout|| = ",ES12.5)') norm 402 end if 403 call DMRestoreGlobalVector(dmUA2, tmpVec,ierr);CHKERRA(ierr) 404 405 ! Building and saving Sigma 406 ! We set sigma_0 = rank (to see partitioning) 407 ! sigma_1 = cell set ID 408 ! sigma_2 = x_coordinate of the cell center of mass 409 call DMGetCoordinateSection(dmS, coordSection,ierr);CHKERRA(ierr) 410 call DMGetCoordinatesLocal(dmS, coord,ierr);CHKERRA(ierr) 411 call DMGetLabelIdIS(dmS, "Cell Sets", csIS,ierr);CHKERRA(ierr) 412 call DMGetLabelSize(dmS, "Cell Sets",numCS,ierr);CHKERRA(ierr) 413 call ISGetIndicesF90(csIS, csID,ierr);CHKERRA(ierr) 414 415 do set = 1, numCS 416 call DMGetStratumIS(dmS, "Cell Sets", csID(set), cellIS,ierr);CHKERRA(ierr) 417 call ISGetIndicesF90(cellIS, cellID,ierr);CHKERRA(ierr) 418 call ISGetSize(cellIS, numCells,ierr);CHKERRA(ierr) 419 do cell = 1,numCells 420 call DMPlexVecGetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval,ierr);CHKERRA(ierr) 421 call DMPlexVecGetClosure(dmS, coordSection, coord, cellID(cell), xyz,ierr);CHKERRA(ierr) 422 cval(1) = rank 423 cval(2) = csID(set) 424 cval(3) = 0.0_kPR 425 do c = 1, size(xyz),sdim 426 cval(3) = cval(3) + xyz(c) 427 end do 428 cval(3) = cval(3) * sdim / size(xyz) 429 call DMPlexVecSetClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval, INSERT_ALL_VALUES,ierr);CHKERRA(ierr) 430 call DMPlexVecRestoreClosure(dmS, PETSC_NULL_SECTION, S, cellID(cell), cval,ierr);CHKERRA(ierr) 431 call DMPlexVecRestoreClosure(dmS, coordSection, coord, cellID(cell), xyz,ierr);CHKERRA(ierr) 432 end do 433 call ISRestoreIndicesF90(cellIS, cellID,ierr);CHKERRA(ierr) 434 call ISDestroy(cellIS,ierr);CHKERRA(ierr) 435 end do 436 call ISRestoreIndicesF90(csIS, csID,ierr);CHKERRA(ierr) 437 call ISDestroy(csIS,ierr);CHKERRA(ierr) 438 call VecViewFromOptions(S, PETSC_NULL_OPTIONS, "-s_vec_view",ierr);CHKERRA(ierr) 439 440 ! Writing zonal variables in Exodus file 441 call DMSetOutputSequenceNumber(dmS,0_kPI,time,ierr);CHKERRA(ierr) 442 call VecView(S,viewer,ierr);CHKERRA(ierr) 443 444 ! Reading zonal variables in Exodus file */ 445 call DMGetGlobalVector(dmS, tmpVec,ierr);CHKERRA(ierr) 446 call VecSet(tmpVec, -1000.0_kPR,ierr);CHKERRA(ierr) 447 call PetscObjectSetName(tmpVec, "Sigma",ierr);CHKERRA(ierr) 448 call VecLoad(tmpVec,viewer,ierr);CHKERRA(ierr) 449 call VecAXPY(S, -1.0_kPR, tmpVec,ierr);CHKERRA(ierr) 450 call VecNorm(S, NORM_INFINITY,norm,ierr);CHKERRQ(ierr) 451 if (norm > PETSC_SQRT_MACHINE_EPSILON) then 452 write(IOBuffer,'("Sigma ||Vin - Vout|| = ",ES12.5)') norm 453 end if 454 call DMRestoreGlobalVector(dmS, tmpVec,ierr);CHKERRA(ierr) 455 456 call DMRestoreGlobalVector(dmUA2, UA2,ierr);CHKERRA(ierr) 457 call DMRestoreGlobalVector(dmUA, UA,ierr);CHKERRA(ierr) 458 call DMRestoreGlobalVector(dmS, S,ierr);CHKERRA(ierr) 459 call DMRestoreGlobalVector(dmA, A,ierr);CHKERRA(ierr) 460 call DMRestoreGlobalVector(dmU, U,ierr);CHKERRA(ierr) 461 call DMRestoreGlobalVector(dm, X,ierr);CHKERRA(ierr) 462 call DMDestroy(dmU,ierr);CHKERRA(ierr); 463 call ISDestroy(isU,ierr);CHKERRA(ierr) 464 call DMDestroy(dmA,ierr);CHKERRA(ierr); 465 call ISDestroy(isA,ierr);CHKERRA(ierr) 466 call DMDestroy(dmS,ierr);CHKERRA(ierr); 467 call ISDestroy(isS,ierr);CHKERRA(ierr) 468 call DMDestroy(dmUA,ierr);CHKERRA(ierr) 469 call ISDestroy(isUA,ierr);CHKERRA(ierr) 470 call DMDestroy(dmUA2,ierr);CHKERRA(ierr) 471 call DMDestroy(dm,ierr);CHKERRA(ierr) 472 473 deallocate(pStartDepth) 474 deallocate(pEndDepth) 475 476 call PetscViewerDestroy(viewer,ierr);CHKERRA(ierr) 477 call PetscFinalize(ierr) 478end program ex26f90 479 480! /*TEST 481! 482! build: 483! requires: exodusii pnetcdf !complex 484! # 2D seq 485! test: 486! suffix: 0 487! 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 488! #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 489! test: 490! suffix: 1 491! 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 492! 493! test: 494! suffix: 2 495! 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 496! #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 497! test: 498! suffix: 3 499! 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 500! test: 501! suffix: 4 502! 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 503! test: 504! suffix: 5 505! 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 506! # 2D par 507! test: 508! suffix: 6 509! nsize: 2 510! 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 511! #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 512! test: 513! suffix: 7 514! nsize: 2 515! 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 516! test: 517! suffix: 8 518! nsize: 2 519! 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 520! #TODO: bug in call to NetCDF failed to complete invalid type definition in file id 65536 NetCDF: invalid dimension ID or name 521! test: 522! suffix: 9 523! nsize: 2 524! 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 525! test: 526! suffix: 10 527! nsize: 2 528! 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 529! test: 530! # Something is now broken with parallel read/write for wahtever shape H is 531! TODO: broken 532! suffix: 11 533! nsize: 2 534! 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 535 536! #3d seq 537! test: 538! suffix: 12 539! 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 540! test: 541! suffix: 13 542! 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 543! #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 544! test: 545! suffix: 14 546! 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 547! test: 548! suffix: 15 549! 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 550! #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 551! #3d par 552! test: 553! suffix: 16 554! nsize: 2 555! 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 556! test: 557! suffix: 17 558! nsize: 2 559! 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 560! #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 561! test: 562! suffix: 18 563! nsize: 2 564! 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 565! test: 566! suffix: 19 567! nsize: 2 568! 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 569! #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 570! TEST*/ 571