Actual source code: ex22f.F
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/petscda.h
29: DMMG dmmg
30: PetscErrorCode ierr
31: DA da
32: external ComputeRHS,ComputeJacobian
33: Vec x
34: PetscInt i1,i3
36: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
38: i3 = 3
39: i1 = 1
40: call DMMGCreate(PETSC_COMM_WORLD,i3,PETSC_NULL_INTEGER,dmmg,ierr)
41: call DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR, &
42: & i3,i3,i3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,i1,i1, &
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 finclude/petscsys.h
65: #include finclude/petscvec.h
66: #include finclude/petscda.h
68: PetscErrorCode ierr
69: PetscInt mx,my,mz
70: PetscScalar h
71: Vec b
72: DMMG dmmg
73: DA da
75: call DMMGGetDA(dmmg,da,ierr)
76: call DAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz, &
77: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
78: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
79: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
80: & PETSC_NULL_INTEGER,ierr)
81: h = 1.d0/((mx-1)*(my-1)*(mz-1))
83: call VecSet(b,h,ierr)
84: return
85: end
87:
88: subroutine ComputeJacobian(dmmg,JJ,jac,ierr)
89: implicit none
91: #include finclude/petscsys.h
92: #include finclude/petscvec.h
93: #include finclude/petscmat.h
94: #include finclude/petscda.h
96: DMMG dmmg
97: Mat jac,JJ
98: PetscErrorCode ierr
100: DA da
101: PetscInt i,j,k,mx,my,mz,xm
102: PetscInt ym,zm,xs,ys,zs,i1,i7
103: PetscScalar v(7),Hx,Hy,Hz
104: PetscScalar HxHydHz,HyHzdHx
105: PetscScalar HxHzdHy
106: MatStencil row(4),col(4,7)
108: i1 = 1
109: i7 = 7
110: call DMMGGetDA(dmmg,da,ierr)
111: call DAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz, &
112: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
113: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
114: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
115: & PETSC_NULL_INTEGER,ierr)
117: Hx = 1.d0 / (mx-1)
118: Hy = 1.d0 / (my-1)
119: Hz = 1.d0 / (mz-1)
120: HxHydHz = Hx*Hy/Hz
121: HxHzdHy = Hx*Hz/Hy
122: HyHzdHx = Hy*Hz/Hx
123: call DAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)
124:
125: do 10,k=zs,zs+zm-1
126: do 20,j=ys,ys+ym-1
127: do 30,i=xs,xs+xm-1
128: row(MatStencil_i) = i
129: row(MatStencil_j) = j
130: row(MatStencil_k) = k
131: if (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or. &
132: & j.eq.my-1 .or. k.eq.mz-1) then
133: v(1) = 2.d0*(HxHydHz + HxHzdHy + HyHzdHx)
134: call MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES, &
135: & ierr)
136: else
137: v(1) = -HxHydHz
138: col(MatStencil_i,1) = i
139: col(MatStencil_j,1) = j
140: col(MatStencil_k,1) = k-1
141: v(2) = -HxHzdHy
142: col(MatStencil_i,2) = i
143: col(MatStencil_j,2) = j-1
144: col(MatStencil_k,2) = k
145: v(3) = -HyHzdHx
146: col(MatStencil_i,3) = i-1
147: col(MatStencil_j,3) = j
148: col(MatStencil_k,3) = k
149: v(4) = 2.d0*(HxHydHz + HxHzdHy + HyHzdHx)
150: col(MatStencil_i,4) = i
151: col(MatStencil_j,4) = j
152: col(MatStencil_k,4) = k
153: v(5) = -HyHzdHx
154: col(MatStencil_i,5) = i+1
155: col(MatStencil_j,5) = j
156: col(MatStencil_k,5) = k
157: v(6) = -HxHzdHy
158: col(MatStencil_i,6) = i
159: col(MatStencil_j,6) = j+1
160: col(MatStencil_k,6) = k
161: v(7) = -HxHydHz
162: col(MatStencil_i,7) = i
163: col(MatStencil_j,7) = j
164: col(MatStencil_k,7) = k+1
165: call MatSetValuesStencil(jac,i1,row,i7,col,v,INSERT_VALUES, &
166: & ierr)
167: endif
168: 30 continue
169: 20 continue
170: 10 continue
172: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
173: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
174: return
175: end