Actual source code: ex62f90.F90
1: program ex62f90
2: #include "petsc/finclude/petsc.h"
3: use petsc
4: implicit none
5: #include "exodusII.inc"
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)
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
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
42: type(tIS) :: cellIS
43: PetscInt,dimension(:),pointer :: cellID
44: PetscInt :: numCells, cell, closureSize
45: PetscInt,dimension(:),pointer :: closureA,closure
47: type(tPetscSection) :: sectionUA,coordSection
48: type(tVec) :: UALoc,coord
49: PetscScalar,dimension(:),pointer :: cval,xyz
50: PetscInt :: dofUA,offUA,c
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
73: type(tPetscSF) :: migrationSF, natSF, natPointSF, natPointSFInv
74: PetscPartitioner :: part
76: type(tVec) :: tmpVec
77: PetscReal :: norm
78: PetscReal :: time = 1.234_kPR
80: PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER,ierr))
81: if (ierr /= 0) then
82: print*,'Unable to initialize PETSc'
83: stop
84: endif
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 ')
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