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