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