module PETSc_Tools #include "petsc/finclude/petscdmplex.h" #include "petsc/finclude/petscdmlabel.h" use petscsys use petscdmplex use IO_Treat, only: Outputs implicit none public ! default as public !----------------------------------! ! General PETSc variables ! !----------------------------------! PetscErrorCode :: ierr PetscMPIInt :: myid, nproc, icpu PetscInt :: size !----------------------------------! ! PETSc section and output field ! !----------------------------------! DMLabel, pointer :: nolabel(:) => NULL() Vec :: node ! node looping vector PetscSection :: section PetscInt :: numFields = 3 ! Note: numFields does not mean "number of variables" ! Example: ! numFields=1: put data on vertices ! numFields=2: put data on cells ! numFields=3: put data on faces PetscInt :: numBC = 1 PetscInt, target, dimension(3) :: numComp ! dont fix size. Its okay for Finite-Volume PetscInt, target, dimension(12) :: numDof ! dont fix size. Its okay for Finite-Volume PetscInt, pointer :: pNumComp(:) PetscInt, pointer :: pNumDof(:) PetscInt, target, dimension(1) :: bcField PetscInt, pointer :: pBcField(:) IS, target, dimension(1) :: bcCompIS IS, target, dimension(1) :: bcPointIS IS, pointer :: pBcCompIS(:) IS, pointer :: pBcPointIS(:) PetscViewer :: viewer !----------------------------------! ! DMPlex variables to handle grid ! !----------------------------------! DM :: dm_geom ! dmplex object encapsulates loaded mesh Vec :: coord ! coordinate array from PETSc PetscInt :: depth, dim ! depth and dimension in spatial grid PetscInt :: pst, pend ! start & end index for entire topological points PetscInt :: vst, vend ! start & end index for verticies PetscInt :: est, eend ! start & end index for edges PetscInt :: cst, cend ! start & end index for cells PetscInt :: fst, fend ! start & end index for faces PetscInt :: nnode ! number of vertices PetscInt :: nedge ! number of edges PetscInt :: ncell ! number of cells ! PetscInt :: nface ! number of faces PetscInt :: nebnd ! number of boundary edge (physical) PetscInt :: nebnda ! number of boundary edge (physical & parallel) PetscInt :: nebndp ! number of boundary edge (parallel) PetscInt :: neinn ! number of inner edge (physical & parallel boundaries excluded) DMLabel :: label_physical ! label for physical boundary (exclude parallel boundary) DMLabel :: label_para ! label for parallel boundary DMLabel :: label_all ! label for all kinds of boundary (both physical and parallel) DMLabel :: label_cell ! label for cell shape PetscInt, parameter :: marker_boundary = 100 ! marker that includes both physical and parallel boundaries PetscInt, parameter :: marker_para = 200 ! marker for parallel boundary !---------------------------------------! ! PETSc objects for halo communication ! !---------------------------------------! PetscSF :: sf!, point_sf, section_sf PetscSection :: s type(tVec) :: vloc ! local vector type(tVec) :: vglo ! global vector !------------------------------------------! ! Dummy variables for vector manipulation ! !------------------------------------------! !-----------------------! ! general for 2D and 3D ! !-----------------------! PetscReal, target, dimension(1) :: sol_tmp PetscReal, pointer :: pSol(:) !=> sol_tmp PetscReal, target, dimension(1) :: xx, yy, zz PetscReal, pointer :: pxx(:) !=> xx PetscReal, pointer :: pyy(:) !=> yy PetscReal, pointer :: pzz(:) !=> zz double precision, allocatable :: sol(:) ! temporary solution vector. !-----------------------! ! for 2D ! !-----------------------! PetscInt, target, dimension(4) :: cone_2D ! number 4 covers triangle and quadrilateral cells in 2D PetscInt, pointer :: pCone_2D(:) !=> cone_2D PetscInt, pointer :: pCone_tmp_2D(:) !=> cone_2D PetscInt, target, dimension(10) :: supp_2D ! number 10 denotes upper limit of supporting edges at given vertex in 2D PetscInt, pointer :: pSupp_2D(:) !=> supp_2D PetscReal, target, dimension(2) :: cent_2D, norm_2D ! number 2 means x and y coordinates PetscReal, pointer :: pCent_2D(:) !=> cent_2D PetscReal, pointer :: pNorm_2D(:) !=> norm_2D contains !----------------------------------! subroutine PETSc_Output_Field() ! print out field data into a vtu file use Utils_Parameters implicit none integer :: i type(Outputs) :: Output character(str_l) :: fname ! Create a scalar field, a vector field, and a surface vector field numComp(:) = 1 ! Note: dont know about exact meaning yet pNumComp => numComp do i=1,numFields*(dim+1) numDof(i) = 0 end do ! define first field on vertices numDof(0*(dim+1)+1) = 1 ! number of content in "vertex" field ! if this number is equal to 1, scalar will be printed ! if this number is lager than 1, array will be printed ! define second field on cells numDof(1*(dim+1)+dim+1) = 0 ! number of content in "cell" field ! Note: even without defining cell field, ! Rank field still printed out ! define third field on faces numDof(2*(dim+1)+dim) = 0 ! number of content in "face" field ! Note: face field will not shown in VTK pNumDof => numDof ! Setup boundary conditions numBC = 1 bcField(1) = 0 pBcField => bcField call ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, bcCompIS(1), ierr);CHKERRA(ierr) pBcCompIS => bcCompIS ! Note : should turn on marker label using option (see run file) ! using other labels failed call DMGetStratumIS(dm_geom, "marker", 1, bcPointIS(1),ierr);CHKERRA(ierr) pBcPointIS => bcPointIS ! Create a PetscSection with this data layout ! nolabel(:) => NULL() call DMSetNumFields(dm_geom, numFields,ierr);CHKERRA(ierr) call DMPlexCreateSection(dm_geom,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr) CHKERRA(ierr) call ISDestroy(bcCompIS(1), ierr);CHKERRA(ierr) call ISDestroy(bcPointIS(1), ierr);CHKERRA(ierr) ! Name the Field variables call PetscSectionSetFieldName(section, 0, 'Node', ierr);CHKERRA(ierr) ! value at node call DMSetLocalSection(dm_geom, section, ierr);CHKERRA(ierr) call PetscSectionDestroy(section, ierr);CHKERRA(ierr) ! Create a Vec with this layout of section ! Note : vector 'node' will include entire fields defined above. call DMGetLocalVector(dm_geom, node, ierr);CHKERRA(ierr) ! Note : Local vector size depends on the way section defined call VecGetSize(node, size, ierr);CHKERRA(ierr) do i=0,size-1 ! vertex first and then face or vertex : according to your section definition ! example - 1 : fill vertex vector and face factor separately ! if(i.le.nnode-1) then ! fill node ! ! do whatever you want, and then you fill solution vector as below before print it out ! ! using section layout we defined above ! call VecSetValues(node, 1, i, one, INSERT_VALUES, ierr);CHKERRA(ierr) ! else ! fill face (or cell center value) ! call VecSetValues(node, 1, i, two, INSERT_VALUES, ierr);CHKERRA(ierr) ! endif ! example - 2 : pass arbitrary non-petsc solution vector to petsc data layout pSol = sol(i+1) call VecSetValues(node, 1, i, pSol, INSERT_VALUES, ierr);CHKERRA(ierr) enddo ! create viewer call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr);CHKERRA(ierr) ! set viewer type call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr);CHKERRA(ierr) ! set output file name for viewer call Output%Path fname = Output%Output_Path_PETSc_Field call PetscViewerFileSetName(viewer, fname, ierr);CHKERRA(ierr) !call PetscViewerFileSetName(viewer, 'sol.vtu', ierr);CHKERRA(ierr) ! print dm object to the file call VecView(node, viewer, ierr);CHKERRA(ierr) ! destroy call PetscViewerDestroy(viewer, ierr);CHKERRA(ierr) call DMRestoreLocalVector(dm_geom, node, ierr);CHKERRA(ierr) call DMDestroy(dm_geom, ierr);CHKERRA(ierr) end subroutine PETSc_Output_Field end module PETSc_Tools