Actual source code: ex13f90.F90
1: !
2: !
3: !/*T
4: ! Concepts: KSP^basic sequential example
5: ! Concepts: KSP^Laplacian, 2d
6: ! Concepts: Laplacian, 2d
7: ! Processors: 1
8: !T*/
9: ! -----------------------------------------------------------------------
11: module UserModule
12: #include <petsc/finclude/petscksp.h>
13: use petscksp
14: type User
15: Vec x
16: Vec b
17: Mat A
18: KSP ksp
19: PetscInt m
20: PetscInt n
21: end type User
22: end module
24: program main
25: use UserModule
26: implicit none
28: ! User-defined context that contains all the data structures used
29: ! in the linear solution process.
31: ! Vec x,b /* solution vector, right hand side vector and work vector */
32: ! Mat A /* sparse matrix */
33: ! KSP ksp /* linear solver context */
34: ! int m,n /* grid dimensions */
35: !
36: ! Since we cannot store Scalars and integers in the same context,
37: ! we store the integers/pointers in the user-defined context, and
38: ! the scalar values are carried in the common block.
39: ! The scalar values in this simplistic example could easily
40: ! be recalculated in each routine, where they are needed.
41: !
42: ! Scalar hx2,hy2 /* 1/(m+1)*(m+1) and 1/(n+1)*(n+1) */
44: ! Note: Any user-defined Fortran routines MUST be declared as external.
46: external UserInitializeLinearSolver
47: external UserFinalizeLinearSolver
48: external UserDoLinearSolver
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: ! Variable declarations
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: PetscScalar hx,hy,x,y
55: type(User) userctx
56: PetscErrorCode ierr
57: PetscInt m,n,t,tmax,i,j
58: PetscBool flg
59: PetscMPIInt size
60: PetscReal enorm
61: PetscScalar cnorm
62: PetscScalar,ALLOCATABLE :: userx(:,:)
63: PetscScalar,ALLOCATABLE :: userb(:,:)
64: PetscScalar,ALLOCATABLE :: solution(:,:)
65: PetscScalar,ALLOCATABLE :: rho(:,:)
67: PetscReal hx2,hy2
68: common /param/ hx2,hy2
70: tmax = 2
71: m = 6
72: n = 7
74: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: ! Beginning of program
76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
79: if (ierr .ne. 0) then
80: print*,'Unable to initialize PETSc'
81: stop
82: endif
83: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
84: if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
86: ! The next two lines are for testing only; these allow the user to
87: ! decide the grid size at runtime.
89: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr);CHKERRA(ierr)
90: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr);CHKERRA(ierr)
92: ! Create the empty sparse matrix and linear solver data structures
94: call UserInitializeLinearSolver(m,n,userctx,ierr);CHKERRA(ierr)
96: ! Allocate arrays to hold the solution to the linear system. This
97: ! approach is not normally done in PETSc programs, but in this case,
98: ! since we are calling these routines from a non-PETSc program, we
99: ! would like to reuse the data structures from another code. So in
100: ! the context of a larger application these would be provided by
101: ! other (non-PETSc) parts of the application code.
103: ALLOCATE (userx(m,n),userb(m,n),solution(m,n))
105: ! Allocate an array to hold the coefficients in the elliptic operator
107: ALLOCATE (rho(m,n))
109: ! Fill up the array rho[] with the function rho(x,y) = x; fill the
110: ! right-hand-side b[] and the solution with a known problem for testing.
112: hx = 1.0/real(m+1)
113: hy = 1.0/real(n+1)
114: y = hy
115: do 20 j=1,n
116: x = hx
117: do 10 i=1,m
118: rho(i,j) = x
119: solution(i,j) = sin(2.*PETSC_PI*x)*sin(2.*PETSC_PI*y)
120: userb(i,j) = -2.*PETSC_PI*cos(2.*PETSC_PI*x)*sin(2.*PETSC_PI*y) + &
121: & 8*PETSC_PI*PETSC_PI*x*sin(2.*PETSC_PI*x)*sin(2.*PETSC_PI*y)
122: x = x + hx
123: 10 continue
124: y = y + hy
125: 20 continue
127: ! Loop over a bunch of timesteps, setting up and solver the linear
128: ! system for each time-step.
129: ! Note that this loop is somewhat artificial. It is intended to
130: ! demonstrate how one may reuse the linear solvers in each time-step.
132: do 100 t=1,tmax
133: call UserDoLinearSolver(rho,userctx,userb,userx,ierr);CHKERRA(ierr)
135: ! Compute error: Note that this could (and usually should) all be done
136: ! using the PETSc vector operations. Here we demonstrate using more
137: ! standard programming practices to show how they may be mixed with
138: ! PETSc.
139: cnorm = 0.0
140: do 90 j=1,n
141: do 80 i=1,m
142: cnorm = cnorm + PetscConj(solution(i,j)-userx(i,j))*(solution(i,j)-userx(i,j))
143: 80 continue
144: 90 continue
145: enorm = PetscRealPart(cnorm*hx*hy)
146: write(6,115) m,n,enorm
147: 115 format ('m = ',I2,' n = ',I2,' error norm = ',1PE11.4)
148: 100 continue
150: ! We are finished solving linear systems, so we clean up the
151: ! data structures.
153: DEALLOCATE (userx,userb,solution,rho)
155: call UserFinalizeLinearSolver(userctx,ierr);CHKERRA(ierr)
156: call PetscFinalize(ierr)
157: end
159: ! ----------------------------------------------------------------
160: subroutine UserInitializeLinearSolver(m,n,userctx,ierr)
161: use UserModule
162: implicit none
164: PetscInt m,n
165: PetscErrorCode ierr
166: type(User) userctx
168: common /param/ hx2,hy2
169: PetscReal hx2,hy2
171: ! Local variable declararions
172: Mat A
173: Vec b,x
174: KSP ksp
175: PetscInt Ntot,five,one
177: ! Here we assume use of a grid of size m x n, with all points on the
178: ! interior of the domain, i.e., we do not include the points corresponding
179: ! to homogeneous Dirichlet boundary conditions. We assume that the domain
180: ! is [0,1]x[0,1].
182: hx2 = (m+1)*(m+1)
183: hy2 = (n+1)*(n+1)
184: Ntot = m*n
186: five = 5
187: one = 1
189: ! Create the sparse matrix. Preallocate 5 nonzeros per row.
191: call MatCreateSeqAIJ(PETSC_COMM_SELF,Ntot,Ntot,five,PETSC_NULL_INTEGER,A,ierr);PetscCall(ierr)
192: !
193: ! Create vectors. Here we create vectors with no memory allocated.
194: ! This way, we can use the data structures already in the program
195: ! by using VecPlaceArray() subroutine at a later stage.
196: !
197: call VecCreateSeqWithArray(PETSC_COMM_SELF,one,Ntot,PETSC_NULL_SCALAR,b,ierr);PetscCall(ierr)
198: call VecDuplicate(b,x,ierr);PetscCall(ierr)
200: ! Create linear solver context. This will be used repeatedly for all
201: ! the linear solves needed.
203: call KSPCreate(PETSC_COMM_SELF,ksp,ierr);PetscCall(ierr)
205: userctx%x = x
206: userctx%b = b
207: userctx%A = A
208: userctx%ksp = ksp
209: userctx%m = m
210: userctx%n = n
212: return
213: end
214: ! -----------------------------------------------------------------------
216: ! Solves -div (rho grad psi) = F using finite differences.
217: ! rho is a 2-dimensional array of size m by n, stored in Fortran
218: ! style by columns. userb is a standard one-dimensional array.
220: subroutine UserDoLinearSolver(rho,userctx,userb,userx,ierr)
221: use UserModule
222: implicit none
224: PetscErrorCode ierr
225: type(User) userctx
226: PetscScalar rho(*),userb(*),userx(*)
228: common /param/ hx2,hy2
229: PetscReal hx2,hy2
231: PC pc
232: KSP ksp
233: Vec b,x
234: Mat A
235: PetscInt m,n,one
236: PetscInt i,j,II,JJ
237: PetscScalar v
239: one = 1
240: x = userctx%x
241: b = userctx%b
242: A = userctx%A
243: ksp = userctx%ksp
244: m = userctx%m
245: n = userctx%n
247: ! This is not the most efficient way of generating the matrix,
248: ! but let's not worry about it. We should have separate code for
249: ! the four corners, each edge and then the interior. Then we won't
250: ! have the slow if-tests inside the loop.
251: !
252: ! Compute the operator
253: ! -div rho grad
254: ! on an m by n grid with zero Dirichlet boundary conditions. The rho
255: ! is assumed to be given on the same grid as the finite difference
256: ! stencil is applied. For a staggered grid, one would have to change
257: ! things slightly.
259: II = 0
260: do 110 j=1,n
261: do 100 i=1,m
262: if (j .gt. 1) then
263: JJ = II - m
264: v = -0.5*(rho(II+1) + rho(JJ+1))*hy2
265: call MatSetValues(A,one,II,one,JJ,v,INSERT_VALUES,ierr);PetscCall(ierr)
266: endif
267: if (j .lt. n) then
268: JJ = II + m
269: v = -0.5*(rho(II+1) + rho(JJ+1))*hy2
270: call MatSetValues(A,one,II,one,JJ,v,INSERT_VALUES,ierr);PetscCall(ierr)
271: endif
272: if (i .gt. 1) then
273: JJ = II - 1
274: v = -0.5*(rho(II+1) + rho(JJ+1))*hx2
275: call MatSetValues(A,one,II,one,JJ,v,INSERT_VALUES,ierr);PetscCall(ierr)
276: endif
277: if (i .lt. m) then
278: JJ = II + 1
279: v = -0.5*(rho(II+1) + rho(JJ+1))*hx2
280: call MatSetValues(A,one,II,one,JJ,v,INSERT_VALUES,ierr);PetscCall(ierr)
281: endif
282: v = 2*rho(II+1)*(hx2+hy2)
283: call MatSetValues(A,one,II,one,II,v,INSERT_VALUES,ierr);PetscCall(ierr)
284: II = II+1
285: 100 continue
286: 110 continue
287: !
288: ! Assemble matrix
289: !
290: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);PetscCall(ierr)
291: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);PetscCall(ierr)
293: !
294: ! Set operators. Here the matrix that defines the linear system
295: ! also serves as the preconditioning matrix. Since all the matrices
296: ! will have the same nonzero pattern here, we indicate this so the
297: ! linear solvers can take advantage of this.
298: !
299: call KSPSetOperators(ksp,A,A,ierr);PetscCall(ierr)
301: !
302: ! Set linear solver defaults for this problem (optional).
303: ! - Here we set it to use direct LU factorization for the solution
304: !
305: call KSPGetPC(ksp,pc,ierr);PetscCall(ierr)
306: call PCSetType(pc,PCLU,ierr);PetscCall(ierr)
308: !
309: ! Set runtime options, e.g.,
310: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
311: ! These options will override those specified above as long as
312: ! KSPSetFromOptions() is called _after_ any other customization
313: ! routines.
314: !
315: ! Run the program with the option -help to see all the possible
316: ! linear solver options.
317: !
318: call KSPSetFromOptions(ksp,ierr);PetscCall(ierr)
320: !
321: ! This allows the PETSc linear solvers to compute the solution
322: ! directly in the user's array rather than in the PETSc vector.
323: !
324: ! This is essentially a hack and not highly recommend unless you
325: ! are quite comfortable with using PETSc. In general, users should
326: ! write their entire application using PETSc vectors rather than
327: ! arrays.
328: !
329: call VecPlaceArray(x,userx,ierr);PetscCall(ierr)
330: call VecPlaceArray(b,userb,ierr);PetscCall(ierr)
332: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333: ! Solve the linear system
334: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
336: call KSPSolve(ksp,b,x,ierr);PetscCall(ierr)
338: call VecResetArray(x,ierr);PetscCall(ierr)
339: call VecResetArray(b,ierr);PetscCall(ierr)
341: return
342: end
344: ! ------------------------------------------------------------------------
346: subroutine UserFinalizeLinearSolver(userctx,ierr)
347: use UserModule
348: implicit none
350: PetscErrorCode ierr
351: type(User) userctx
353: !
354: ! We are all done and don't need to solve any more linear systems, so
355: ! we free the work space. All PETSc objects should be destroyed when
356: ! they are no longer needed.
357: !
358: call VecDestroy(userctx%x,ierr);PetscCall(ierr)
359: call VecDestroy(userctx%b,ierr);PetscCall(ierr)
360: call MatDestroy(userctx%A,ierr);PetscCall(ierr)
361: call KSPDestroy(userctx%ksp,ierr);PetscCall(ierr)
363: return
364: end
366: !
367: !/*TEST
368: !
369: ! test:
370: ! args: -m 19 -n 20 -ksp_gmres_cgs_refinement_type refine_always
371: ! output_file: output/ex13f90_1.out
372: !
373: !TEST*/