Actual source code: ex22f.F
petsc-3.3-p6 2013-02-11
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
36: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
38: i3 = -3
39: i1 = 1
40: call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
41: call DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, &
42: & DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE, &
43: & DMDA_STENCIL_STAR,i3,i3,i3,PETSC_DECIDE,PETSC_DECIDE, &
44: & PETSC_DECIDE,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
45: & PETSC_NULL_INTEGER,da,ierr)
46: call DMSetFunction(da,ComputeRHS,ierr)
47: call DMSetJacobian(da,ComputeMatrix,ierr)
48: call KSPSetDM(ksp,da,ierr)
50: call KSPSetFromOptions(ksp,ierr)
51: call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
52: call KSPGetSolution(ksp,x,ierr)
53: ! call VecView(x,PETSC_NULL_OBJECT,ierr)
54: call KSPDestroy(ksp,ierr)
55: call DMDestroy(da,ierr)
56: call PetscFinalize(ierr)
58: end
60: subroutine ComputeRHS(da,x,b,ierr)
61: implicit none
63: #include <finclude/petscsys.h>
64: #include <finclude/petscvec.h>
65: #include <finclude/petscdmda.h>
67: PetscErrorCode ierr
68: PetscInt mx,my,mz
69: PetscScalar h
70: Vec x,b
71: DM da
73: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz, &
74: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
75: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
76: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
77: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
78: & PETSC_NULL_INTEGER,ierr)
79: h = 1.d0/((mx-1)*(my-1)*(mz-1))
81: call VecSet(b,h,ierr)
82: return
83: end
85:
86: subroutine ComputeMatrix(da,x,JJ,jac,str,ierr)
87: implicit none
89: #include <finclude/petscsys.h>
90: #include <finclude/petscvec.h>
91: #include <finclude/petscmat.h>
92: #include <finclude/petscdmda.h>
94: Mat jac,JJ
95: PetscErrorCode ierr
97: DM da
98: PetscInt i,j,k,mx,my,mz,xm
99: PetscInt ym,zm,xs,ys,zs,i1,i7
100: PetscScalar v(7),Hx,Hy,Hz
101: PetscScalar HxHydHz,HyHzdHx
102: PetscScalar HxHzdHy
103: MatStencil row(4),col(4,7)
104: Vec x
105: MatStructure str
106: i1 = 1
107: i7 = 7
108: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz, &
109: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
110: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
111: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
112: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
113: & PETSC_NULL_INTEGER,ierr)
115: Hx = 1.d0 / (mx-1)
116: Hy = 1.d0 / (my-1)
117: Hz = 1.d0 / (mz-1)
118: HxHydHz = Hx*Hy/Hz
119: HxHzdHy = Hx*Hz/Hy
120: HyHzdHx = Hy*Hz/Hx
121: call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)
122:
123: do 10,k=zs,zs+zm-1
124: do 20,j=ys,ys+ym-1
125: do 30,i=xs,xs+xm-1
126: row(MatStencil_i) = i
127: row(MatStencil_j) = j
128: row(MatStencil_k) = k
129: if (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or. &
130: & j.eq.my-1 .or. k.eq.mz-1) then
131: v(1) = 2.d0*(HxHydHz + HxHzdHy + HyHzdHx)
132: call MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES, &
133: & ierr)
134: else
135: v(1) = -HxHydHz
136: col(MatStencil_i,1) = i
137: col(MatStencil_j,1) = j
138: col(MatStencil_k,1) = k-1
139: v(2) = -HxHzdHy
140: col(MatStencil_i,2) = i
141: col(MatStencil_j,2) = j-1
142: col(MatStencil_k,2) = k
143: v(3) = -HyHzdHx
144: col(MatStencil_i,3) = i-1
145: col(MatStencil_j,3) = j
146: col(MatStencil_k,3) = k
147: v(4) = 2.d0*(HxHydHz + HxHzdHy + HyHzdHx)
148: col(MatStencil_i,4) = i
149: col(MatStencil_j,4) = j
150: col(MatStencil_k,4) = k
151: v(5) = -HyHzdHx
152: col(MatStencil_i,5) = i+1
153: col(MatStencil_j,5) = j
154: col(MatStencil_k,5) = k
155: v(6) = -HxHzdHy
156: col(MatStencil_i,6) = i
157: col(MatStencil_j,6) = j+1
158: col(MatStencil_k,6) = k
159: v(7) = -HxHydHz
160: col(MatStencil_i,7) = i
161: col(MatStencil_j,7) = j
162: col(MatStencil_k,7) = k+1
163: call MatSetValuesStencil(jac,i1,row,i7,col,v,INSERT_VALUES, &
164: & ierr)
165: endif
166: 30 continue
167: 20 continue
168: 10 continue
170: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
171: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
172: return
173: end