Actual source code: ex22f.F

petsc-3.6.3 2015-12-03
Report Typos and Errors
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 <petsc/finclude/petscsys.h>
24: #include <petsc/finclude/petscvec.h>
25: #include <petsc/finclude/petscmat.h>
26: #include <petsc/finclude/petscpc.h>
27: #include <petsc/finclude/petscksp.h>
28: #include <petsc/finclude/petscdm.h>
29: #include <petsc/finclude/petscdmda.h>
30:       PetscErrorCode   ierr
31:       DM               da
32:       KSP              ksp
33:       Vec              x
34:       external         ComputeRHS,ComputeMatrix
35:       PetscInt i1,i3
36:       PetscInt         ctx

38:       call  PetscInitialize(PETSC_NULL_CHARACTER,ierr)

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

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

60:       end

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

65: #include <petsc/finclude/petscsys.h>
66: #include <petsc/finclude/petscvec.h>
67: #include <petsc/finclude/petscdm.h>
68: #include <petsc/finclude/petscdmda.h>

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

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

87:       call VecSet(b,h,ierr)
88:       return
89:       end

92:       subroutine ComputeMatrix(ksp,JJ,jac,ctx,ierr)
93:       implicit none

95: #include <petsc/finclude/petscsys.h>
96: #include <petsc/finclude/petscvec.h>
97: #include <petsc/finclude/petscmat.h>
98: #include <petsc/finclude/petscdm.h>
99: #include <petsc/finclude/petscdmda.h>

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

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

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

177:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
178:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
179:       return
180:       end