Actual source code: ex22f.F
petsc-3.7.5 2017-01-01
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.0/real((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.0 / real(mx-1)
123: Hy = 1.0 / real(my-1)
124: Hz = 1.0 / real(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.0*(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.0*(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