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