Actual source code: ex22f.F

  1: ! $Id: ex22f.F,v 1.6 2001/09/11 15:09:12 bsmith Exp $
  2: !
  3: !   Laplacian in 3D. Modeled by the partial differential equation
  4: !
  5: !   Laplacian u = 1,0 < x,y,z < 1,
  6: !
  7: !   with boundary conditions
  8: !
  9: !   u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
 10: !
 11: !   This uses multigrid to solve the linear system

 13:       program main
 14:       implicit none

 16: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 17: !                    Include files
 18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 19: !
 20: !     petsc.h  - base PETSc routines      petscvec.h - vectors
 21: !     petscsys.h    - system routines          petscmat.h - matrices
 22: !     petscksp.h    - Krylov subspace methods  petscpc.h  - preconditioners
 23: !     petscsles.h   - SLES interface

 25:  #include include/finclude/petsc.h
 26:  #include include/finclude/petscvec.h
 27:  #include include/finclude/petscmat.h
 28:  #include include/finclude/petscpc.h
 29:  #include include/finclude/petscksp.h
 30:  #include include/finclude/petscda.h
 31:       DMMG             dmmg
 32:       integer          ierr
 33:       double precision norm
 34:       DA               da
 35:       external         ComputeRHS,ComputeJacobian
 36:       Vec              x

 38:       call  PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 40:       call DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL_INTEGER,dmmg,ierr)
 41:       call DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,3,            &
 42:      &                3,3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,               &
 43:      &                PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                        &
 44:      &                PETSC_NULL_INTEGER,da,ierr)
 45:       call DMMGSetDM(dmmg,da,ierr)
 46:       call DADestroy(da,ierr)


 49:       call DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian,ierr)

 51:       call DMMGSolve(dmmg,ierr)

 53:       call DMMGGetx(dmmg,x,ierr)
 54: !      call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)

 56:       call DMMGDestroy(dmmg,ierr)
 57:       call PetscFinalize(ierr)

 59:       end

 61:       subroutine ComputeRHS(dmmg,b,ierr)
 62:       implicit none

 64:  #include include/finclude/petsc.h
 65:  #include include/finclude/petscvec.h
 66:  #include include/finclude/petscda.h

 68:       integer      ierr,mx,my,mz
 69:       PetscScalar  h
 70:       Vec          b
 71:       DMMG         dmmg
 72:       DA           da

 74:       call DMMGGetDA(dmmg,da,ierr)
 75:       call DAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,                        &
 76:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 77:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 78:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 79:      &               PETSC_NULL_INTEGER,ierr)
 80:       h    = 1.d0/((mx-1)*(my-1)*(mz-1))

 82:       call VecSet(h,b,ierr)
 83:       return
 84:       end

 86: 
 87:       subroutine ComputeJacobian(dmmg,jac,ierr)
 88:       implicit none

 90:  #include include/finclude/petsc.h
 91:  #include include/finclude/petscvec.h
 92:  #include include/finclude/petscmat.h
 93:  #include include/finclude/petscda.h

 95:       DMMG         dmmg
 96:       Mat          jac
 97:       integer      ierr

 99:       DA           da
100:       integer      i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs
101:       PetscScalar  v(7),Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy
102:       MatStencil   row(4),col(4,7)

104:       call DMMGGetDA(dmmg,da,ierr)
105:       call DAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,                       &
106:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
107:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
108:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
109:      &               PETSC_NULL_INTEGER,ierr)

111:       Hx = 1.d0 / (mx-1)
112:       Hy = 1.d0 / (my-1)
113:       Hz = 1.d0 / (mz-1)
114:       HxHydHz = Hx*Hy/Hz
115:       HxHzdHy = Hx*Hz/Hy
116:       HyHzdHx = Hy*Hz/Hx
117:       call DAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)
118: 
119:       do 10,k=zs,zs+zm-1
120:         do 20,j=ys,ys+ym-1
121:           do 30,i=xs,xs+xm-1
122:           row(MatStencil_i) = i
123:           row(MatStencil_j) = j
124:           row(MatStencil_k) = k
125:           if (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or.         &
126:      &         j.eq.my-1 .or. k.eq.mz-1) then
127:             v(1) = 2.d0*(HxHydHz + HxHzdHy + HyHzdHx)
128:             call MatSetValuesStencil(jac,1,row,1,row,v,INSERT_VALUES,    &
129:      &                               ierr)
130:           else
131:             v(1) = -HxHydHz
132:              col(MatStencil_i,1) = i
133:              col(MatStencil_j,1) = j
134:              col(MatStencil_k,1) = k-1
135:             v(2) = -HxHzdHy
136:              col(MatStencil_i,2) = i
137:              col(MatStencil_j,2) = j-1
138:              col(MatStencil_k,2) = k
139:             v(3) = -HyHzdHx
140:              col(MatStencil_i,3) = i-1
141:              col(MatStencil_j,3) = j
142:              col(MatStencil_k,3) = k
143:             v(4) = 2.d0*(HxHydHz + HxHzdHy + HyHzdHx)
144:              col(MatStencil_i,4) = i
145:              col(MatStencil_j,4) = j
146:              col(MatStencil_k,4) = k
147:             v(5) = -HyHzdHx
148:              col(MatStencil_i,5) = i+1
149:              col(MatStencil_j,5) = j
150:              col(MatStencil_k,5) = k
151:             v(6) = -HxHzdHy
152:              col(MatStencil_i,6) = i
153:              col(MatStencil_j,6) = j+1
154:              col(MatStencil_k,6) = k
155:             v(7) = -HxHydHz
156:              col(MatStencil_i,7) = i
157:              col(MatStencil_j,7) = j
158:              col(MatStencil_k,7) = k+1
159:             call MatSetValuesStencil(jac,1,row,7,col,v,INSERT_VALUES,ierr)
160:           endif
161:  30       continue
162:  20     continue
163:  10   continue

165:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
166:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
167:       return
168:       end