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