Actual source code: ex22f.F90
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
14: #include <petsc/finclude/petscksp.h>
15: use petscdmda
16: use petscksp
17: implicit none
19: PetscErrorCode ierr
20: DM da
21: KSP ksp
22: Vec x
23: external ComputeRHS,ComputeMatrix
24: PetscInt i1,i3
25: PetscInt ctx
27: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
28: if (ierr .ne. 0) then
29: print*,'Unable to initialize PETSc'
30: stop
31: endif
33: i3 = 3
34: i1 = 1
35: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
36: call DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE, &
37: & DMDA_STENCIL_STAR,i3,i3,i3,PETSC_DECIDE,PETSC_DECIDE, &
38: & PETSC_DECIDE,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
39: call DMSetFromOptions(da,ierr)
40: call DMSetUp(da,ierr)
41: call KSPSetDM(ksp,da,ierr)
42: call KSPSetComputeRHS(ksp,ComputeRHS,ctx,ierr)
43: call KSPSetComputeOperators(ksp,ComputeMatrix,ctx,ierr)
45: call KSPSetFromOptions(ksp,ierr)
46: call KSPSolve(ksp,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
47: call KSPGetSolution(ksp,x,ierr)
48: call KSPDestroy(ksp,ierr)
49: call DMDestroy(da,ierr)
50: call PetscFinalize(ierr)
52: end
54: subroutine ComputeRHS(ksp,b,ctx,ierr)
55: use petscksp
56: implicit none
58: PetscErrorCode ierr
59: PetscInt mx,my,mz
60: PetscScalar h
61: Vec b
62: KSP ksp
63: DM da
64: PetscInt ctx
66: call KSPGetDM(ksp,da,ierr)
67: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
68: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
69: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
70: h = 1.0/real((mx-1)*(my-1)*(mz-1))
72: call VecSet(b,h,ierr)
73: return
74: end
76: subroutine ComputeMatrix(ksp,JJ,jac,ctx,ierr)
77: use petscksp
78: implicit none
80: Mat jac,JJ
81: PetscErrorCode ierr
82: KSP ksp
83: DM da
84: PetscInt i,j,k,mx,my,mz,xm
85: PetscInt ym,zm,xs,ys,zs,i1,i7
86: PetscScalar v(7),Hx,Hy,Hz
87: PetscScalar HxHydHz,HyHzdHx
88: PetscScalar HxHzdHy
89: MatStencil row(4),col(4,7)
90: PetscInt ctx
91: i1 = 1
92: i7 = 7
93: call KSPGetDM(ksp,da,ierr)
94: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
95: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
96: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
98: Hx = 1.0 / real(mx-1)
99: Hy = 1.0 / real(my-1)
100: Hz = 1.0 / real(mz-1)
101: HxHydHz = Hx*Hy/Hz
102: HxHzdHy = Hx*Hz/Hy
103: HyHzdHx = Hy*Hz/Hx
104: call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)
106: do 10,k=zs,zs+zm-1
107: do 20,j=ys,ys+ym-1
108: do 30,i=xs,xs+xm-1
109: row(MatStencil_i) = i
110: row(MatStencil_j) = j
111: row(MatStencil_k) = k
112: if (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1 .or. k.eq.mz-1) then
113: v(1) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
114: call MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES,ierr)
115: else
116: v(1) = -HxHydHz
117: col(MatStencil_i,1) = i
118: col(MatStencil_j,1) = j
119: col(MatStencil_k,1) = k-1
120: v(2) = -HxHzdHy
121: col(MatStencil_i,2) = i
122: col(MatStencil_j,2) = j-1
123: col(MatStencil_k,2) = k
124: v(3) = -HyHzdHx
125: col(MatStencil_i,3) = i-1
126: col(MatStencil_j,3) = j
127: col(MatStencil_k,3) = k
128: v(4) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
129: col(MatStencil_i,4) = i
130: col(MatStencil_j,4) = j
131: col(MatStencil_k,4) = k
132: v(5) = -HyHzdHx
133: col(MatStencil_i,5) = i+1
134: col(MatStencil_j,5) = j
135: col(MatStencil_k,5) = k
136: v(6) = -HxHzdHy
137: col(MatStencil_i,6) = i
138: col(MatStencil_j,6) = j+1
139: col(MatStencil_k,6) = k
140: v(7) = -HxHydHz
141: col(MatStencil_i,7) = i
142: col(MatStencil_j,7) = j
143: col(MatStencil_k,7) = k+1
144: call MatSetValuesStencil(jac,i1,row,i7,col,v,INSERT_VALUES,ierr)
145: endif
146: 30 continue
147: 20 continue
148: 10 continue
150: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
151: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
152: return
153: end
155: !/*TEST
156: !
157: ! test:
158: ! args: -pc_mg_type full -ksp_monitor_short -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type preconditioned -pc_type mg -da_refine 2 -ksp_type fgmres
159: ! requires: !single
160: ! output_file: output/ex22_1.out
161: !
162: !TEST*/