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