Actual source code: ex22f.F

petsc-3.4.4 2014-03-13
  1: !
  2: !   Laplacian in 3D. Modeled by the partial differential equation
  3: !
  4: !   Laplacian u = 1,0 < x,y,z < 1,
  5: !
  6: !   with boundary conditions
  7: !
  8: !   u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
  9: !
 10: !   This uses multigrid to solve the linear system

 12:       program main
 13:       implicit none

 15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 16: !                    Include files
 17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 18: !
 19: !     petscsys.h  - base PETSc routines      petscvec.h - vectors
 20: !     petscmat.h - matrices
 21: !     petscksp.h    - Krylov subspace methods  petscpc.h  - preconditioners

 23: #include <finclude/petscsys.h>
 24: #include <finclude/petscvec.h>
 25: #include <finclude/petscmat.h>
 26: #include <finclude/petscpc.h>
 27: #include <finclude/petscksp.h>
 28: #include <finclude/petscdmda.h>
 29:       PetscErrorCode   ierr
 30:       DM               da
 31:       KSP              ksp
 32:       Vec              x
 33:       external         ComputeRHS,ComputeMatrix
 34:       PetscInt i1,i3
 35:       PetscInt         ctx

 37:       call  PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 39:       i3 = -3
 40:       i1 = 1
 41:       call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
 42:       call DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,            &
 43:      &    DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,                               &
 44:      &    DMDA_STENCIL_STAR,i3,i3,i3,PETSC_DECIDE,PETSC_DECIDE,                        &
 45:      &    PETSC_DECIDE,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                        &
 46:      &                PETSC_NULL_INTEGER,da,ierr)
 47:       call KSPSetDM(ksp,da,ierr)
 48:       call KSPSetComputeRHS(ksp,ComputeRHS,ctx,ierr)
 49:       call KSPSetComputeOperators(ksp,ComputeMatrix,ctx,ierr)

 51:       call KSPSetFromOptions(ksp,ierr)
 52:       call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
 53:       call KSPGetSolution(ksp,x,ierr)
 54: !      call VecView(x,PETSC_NULL_OBJECT,ierr)
 55:       call KSPDestroy(ksp,ierr)
 56:       call DMDestroy(da,ierr)
 57:       call PetscFinalize(ierr)

 59:       end

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

 64: #include <finclude/petscsys.h>
 65: #include <finclude/petscvec.h>
 66: #include <finclude/petscdmda.h>

 68:       PetscErrorCode  ierr
 69:       PetscInt mx,my,mz
 70:       PetscScalar  h
 71:       Vec          b
 72:       KSP          ksp
 73:       DM           da
 74:       PetscInt     ctx

 76:       call KSPGetDM(ksp,da,ierr)
 77:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,                        &
 78:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 79:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 80:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 81:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                 &
 82:      &               PETSC_NULL_INTEGER,ierr)
 83:       h    = 1.d0/((mx-1)*(my-1)*(mz-1))

 85:       call VecSet(b,h,ierr)
 86:       return
 87:       end


 90:       subroutine ComputeMatrix(ksp,JJ,jac,str,ctx,ierr)
 91:       implicit none

 93: #include <finclude/petscsys.h>
 94: #include <finclude/petscvec.h>
 95: #include <finclude/petscmat.h>
 96: #include <finclude/petscdmda.h>

 98:       Mat          jac,JJ
 99:       PetscErrorCode    ierr
100:       KSP          ksp
101:       DM           da
102:       PetscInt    i,j,k,mx,my,mz,xm
103:       PetscInt    ym,zm,xs,ys,zs,i1,i7
104:       PetscScalar  v(7),Hx,Hy,Hz
105:       PetscScalar  HxHydHz,HyHzdHx
106:       PetscScalar  HxHzdHy
107:       MatStencil   row(4),col(4,7)
108:       MatStructure str
109:       PetscInt     ctx
110:       i1 = 1
111:       i7 = 7
112:       call KSPGetDM(ksp,da,ierr)
113:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,                       &
114:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
115:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
116:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
117:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                &
118:      &               PETSC_NULL_INTEGER,ierr)

120:       Hx = 1.d0 / (mx-1)
121:       Hy = 1.d0 / (my-1)
122:       Hz = 1.d0 / (mz-1)
123:       HxHydHz = Hx*Hy/Hz
124:       HxHzdHy = Hx*Hz/Hy
125:       HyHzdHx = Hy*Hz/Hx
126:       call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)

128:       do 10,k=zs,zs+zm-1
129:         do 20,j=ys,ys+ym-1
130:           do 30,i=xs,xs+xm-1
131:           row(MatStencil_i) = i
132:           row(MatStencil_j) = j
133:           row(MatStencil_k) = k
134:           if (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or.         &
135:      &         j.eq.my-1 .or. k.eq.mz-1) then
136:             v(1) = 2.d0*(HxHydHz + HxHzdHy + HyHzdHx)
137:             call MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES,    &
138:      &                               ierr)
139:           else
140:             v(1) = -HxHydHz
141:              col(MatStencil_i,1) = i
142:              col(MatStencil_j,1) = j
143:              col(MatStencil_k,1) = k-1
144:             v(2) = -HxHzdHy
145:              col(MatStencil_i,2) = i
146:              col(MatStencil_j,2) = j-1
147:              col(MatStencil_k,2) = k
148:             v(3) = -HyHzdHx
149:              col(MatStencil_i,3) = i-1
150:              col(MatStencil_j,3) = j
151:              col(MatStencil_k,3) = k
152:             v(4) = 2.d0*(HxHydHz + HxHzdHy + HyHzdHx)
153:              col(MatStencil_i,4) = i
154:              col(MatStencil_j,4) = j
155:              col(MatStencil_k,4) = k
156:             v(5) = -HyHzdHx
157:              col(MatStencil_i,5) = i+1
158:              col(MatStencil_j,5) = j
159:              col(MatStencil_k,5) = k
160:             v(6) = -HxHzdHy
161:              col(MatStencil_i,6) = i
162:              col(MatStencil_j,6) = j+1
163:              col(MatStencil_k,6) = k
164:             v(7) = -HxHydHz
165:              col(MatStencil_i,7) = i
166:              col(MatStencil_j,7) = j
167:              col(MatStencil_k,7) = k+1
168:       call MatSetValuesStencil(jac,i1,row,i7,col,v,INSERT_VALUES,               &
169:      &                               ierr)
170:           endif
171:  30       continue
172:  20     continue
173:  10   continue

175:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
176:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
177:       return
178:       end