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