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: program main
12: implicit none
14: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
15: ! Include files
16: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
17: !
18: ! The following include statements are required for KSP Fortran programs:
19: ! petscsys.h - base PETSc routines
20: ! petscvec.h - vectors
21: ! petscmat.h - matrices
22: ! petscksp.h - Krylov subspace methods
23: ! petscpc.h - preconditioners
24: !
25: #include <petsc/finclude/petscsys.h>
26: #include <petsc/finclude/petscvec.h>
27: #include <petsc/finclude/petscmat.h>
28: #include <petsc/finclude/petscksp.h>
29: #include <petsc/finclude/petscpc.h>
31: ! User-defined context that contains all the data structures used
32: ! in the linear solution process.
34: ! Vec x,b /* solution vector, right hand side vector and work vector */
35: ! Mat A /* sparse matrix */
36: ! KSP ksp /* linear solver context */
37: ! int m,n /* grid dimensions */
38: !
39: ! Since we cannot store Scalars and integers in the same context,
40: ! we store the integers/pointers in the user-defined context, and
41: ! the scalar values are carried in the common block.
42: ! The scalar values in this simplistic example could easily
43: ! be recalculated in each routine, where they are needed.
44: !
45: ! Scalar hx2,hy2 /* 1/(m+1)*(m+1) and 1/(n+1)*(n+1) */
47: ! Note: Any user-defined Fortran routines MUST be declared as external.
49: external UserInitializeLinearSolver
50: external UserFinalizeLinearSolver
51: external UserDoLinearSolver
53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: ! Variable declarations
55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: PetscScalar hx,hy,x,y
58: PetscFortranAddr userctx(6)
59: PetscErrorCode ierr
60: PetscInt m,n,t,tmax,i,j
61: PetscBool flg
62: PetscMPIInt size,rank
63: PetscReal enorm
64: PetscScalar cnorm
65: PetscScalar,ALLOCATABLE :: userx(:,:)
66: PetscScalar,ALLOCATABLE :: userb(:,:)
67: PetscScalar,ALLOCATABLE :: solution(:,:)
68: PetscScalar,ALLOCATABLE :: rho(:,:)
70: PetscReal hx2,hy2
71: common /param/ hx2,hy2
73: tmax = 2
74: m = 6
75: n = 7
77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: ! Beginning of program
79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
82: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
83: if (size .ne. 1) then
84: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
85: if (rank .eq. 0) then
86: write(6,*) 'This is a uniprocessor example only!'
87: endif
88: SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)
89: endif
91: ! The next two lines are for testing only; these allow the user to
92: ! decide the grid size at runtime.
94: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
95: & '-m',m,flg,ierr)
96: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
97: & '-n',n,flg,ierr)
99: ! Create the empty sparse matrix and linear solver data structures
101: call UserInitializeLinearSolver(m,n,userctx,ierr)
103: ! Allocate arrays to hold the solution to the linear system. This
104: ! approach is not normally done in PETSc programs, but in this case,
105: ! since we are calling these routines from a non-PETSc program, we
106: ! would like to reuse the data structures from another code. So in
107: ! the context of a larger application these would be provided by
108: ! other (non-PETSc) parts of the application code.
110: ALLOCATE (userx(m,n),userb(m,n),solution(m,n))
112: ! Allocate an array to hold the coefficients in the elliptic operator
114: ALLOCATE (rho(m,n))
116: ! Fill up the array rho[] with the function rho(x,y) = x; fill the
117: ! right-hand-side b[] and the solution with a known problem for testing.
119: hx = 1.0/real(m+1)
120: hy = 1.0/real(n+1)
121: y = hy
122: do 20 j=1,n
123: x = hx
124: do 10 i=1,m
125: rho(i,j) = x
126: solution(i,j) = sin(2.*PETSC_PI*x)*sin(2.*PETSC_PI*y)
127: userb(i,j) = -2.*PETSC_PI*cos(2.*PETSC_PI*x)* &
128: & sin(2.*PETSC_PI*y) + &
129: & 8*PETSC_PI*PETSC_PI*x* &
130: & sin(2.*PETSC_PI*x)*sin(2.*PETSC_PI*y)
131: x = x + hx
132: 10 continue
133: y = y + hy
134: 20 continue
136: ! Loop over a bunch of timesteps, setting up and solver the linear
137: ! system for each time-step.
138: ! Note that this loop is somewhat artificial. It is intended to
139: ! demonstrate how one may reuse the linear solvers in each time-step.
141: do 100 t=1,tmax
142: call UserDoLinearSolver(rho,userctx,userb,userx,ierr)
144: ! Compute error: Note that this could (and usually should) all be done
145: ! using the PETSc vector operations. Here we demonstrate using more
146: ! standard programming practices to show how they may be mixed with
147: ! PETSc.
148: cnorm = 0.0
149: do 90 j=1,n
150: do 80 i=1,m
151: cnorm = cnorm + &
152: & PetscConj(solution(i,j)-userx(i,j))*(solution(i,j)-userx(i,j))
153: 80 continue
154: 90 continue
155: enorm = PetscRealPart(cnorm*hx*hy)
156: write(6,115) m,n,enorm
157: 115 format ('m = ',I2,' n = ',I2,' error norm = ',1PE11.4)
158: 100 continue
160: ! We are finished solving linear systems, so we clean up the
161: ! data structures.
163: DEALLOCATE (userx,userb,solution,rho)
165: call UserFinalizeLinearSolver(userctx,ierr)
166: call PetscFinalize(ierr)
167: end
169: ! ----------------------------------------------------------------
170: subroutine UserInitializeLinearSolver(m,n,userctx,ierr)
172: implicit none
174: #include <petsc/finclude/petscsys.h>
175: #include <petsc/finclude/petscvec.h>
176: #include <petsc/finclude/petscmat.h>
177: #include <petsc/finclude/petscksp.h>
178: #include <petsc/finclude/petscpc.h>
180: PetscInt m,n
181: PetscErrorCode ierr
182: PetscFortranAddr userctx(*)
184: common /param/ hx2,hy2
185: PetscReal hx2,hy2
187: ! Local variable declararions
188: Mat A
189: Vec b,x
190: KSP ksp
191: PetscInt Ntot,five,one
194: ! Here we assume use of a grid of size m x n, with all points on the
195: ! interior of the domain, i.e., we do not include the points corresponding
196: ! to homogeneous Dirichlet boundary conditions. We assume that the domain
197: ! is [0,1]x[0,1].
199: hx2 = (m+1)*(m+1)
200: hy2 = (n+1)*(n+1)
201: Ntot = m*n
203: five = 5
204: one = 1
206: ! Create the sparse matrix. Preallocate 5 nonzeros per row.
208: call MatCreateSeqAIJ(PETSC_COMM_SELF,Ntot,Ntot,five, &
209: & PETSC_NULL_INTEGER,A,ierr)
210: !
211: ! Create vectors. Here we create vectors with no memory allocated.
212: ! This way, we can use the data structures already in the program
213: ! by using VecPlaceArray() subroutine at a later stage.
214: !
215: call VecCreateSeqWithArray(PETSC_COMM_SELF,one,Ntot, &
216: & PETSC_NULL_SCALAR,b,ierr)
217: call VecDuplicate(b,x,ierr)
219: ! Create linear solver context. This will be used repeatedly for all
220: ! the linear solves needed.
222: call KSPCreate(PETSC_COMM_SELF,ksp,ierr)
224: userctx(1) = x
225: userctx(2) = b
226: userctx(3) = A
227: userctx(4) = ksp
228: userctx(5) = m
229: userctx(6) = n
231: return
232: end
233: ! -----------------------------------------------------------------------
235: ! Solves -div (rho grad psi) = F using finite differences.
236: ! rho is a 2-dimensional array of size m by n, stored in Fortran
237: ! style by columns. userb is a standard one-dimensional array.
239: subroutine UserDoLinearSolver(rho,userctx,userb,userx,ierr)
241: implicit none
243: #include <petsc/finclude/petscsys.h>
244: #include <petsc/finclude/petscvec.h>
245: #include <petsc/finclude/petscmat.h>
246: #include <petsc/finclude/petscksp.h>
247: #include <petsc/finclude/petscpc.h>
249: PetscErrorCode ierr
250: PetscFortranAddr userctx(*)
251: PetscScalar rho(*),userb(*),userx(*)
254: common /param/ hx2,hy2
255: PetscReal hx2,hy2
257: PC pc
258: KSP ksp
259: Vec b,x
260: Mat A
261: PetscInt m,n,one
262: PetscInt i,j,II,JJ
263: PetscScalar v
265: ! PetscScalar tmpx(*),tmpb(*)
267: one = 1
268: x = userctx(1)
269: b = userctx(2)
270: A = userctx(3)
271: ksp = userctx(4)
272: m = int(userctx(5))
273: n = int(userctx(6))
275: ! This is not the most efficient way of generating the matrix,
276: ! but let's not worry about it. We should have separate code for
277: ! the four corners, each edge and then the interior. Then we won't
278: ! have the slow if-tests inside the loop.
279: !
280: ! Compute the operator
281: ! -div rho grad
282: ! on an m by n grid with zero Dirichlet boundary conditions. The rho
283: ! is assumed to be given on the same grid as the finite difference
284: ! stencil is applied. For a staggered grid, one would have to change
285: ! things slightly.
287: II = 0
288: do 110 j=1,n
289: do 100 i=1,m
290: if (j .gt. 1) then
291: JJ = II - m
292: v = -0.5*(rho(II+1) + rho(JJ+1))*hy2
293: call MatSetValues(A,one,II,one,JJ,v,INSERT_VALUES,ierr)
294: endif
295: if (j .lt. n) then
296: JJ = II + m
297: v = -0.5*(rho(II+1) + rho(JJ+1))*hy2
298: call MatSetValues(A,one,II,one,JJ,v,INSERT_VALUES,ierr)
299: endif
300: if (i .gt. 1) then
301: JJ = II - 1
302: v = -0.5*(rho(II+1) + rho(JJ+1))*hx2
303: call MatSetValues(A,one,II,one,JJ,v,INSERT_VALUES,ierr)
304: endif
305: if (i .lt. m) then
306: JJ = II + 1
307: v = -0.5*(rho(II+1) + rho(JJ+1))*hx2
308: call MatSetValues(A,one,II,one,JJ,v,INSERT_VALUES,ierr)
309: endif
310: v = 2*rho(II+1)*(hx2+hy2)
311: call MatSetValues(A,one,II,one,II,v,INSERT_VALUES,ierr)
312: II = II+1
313: 100 continue
314: 110 continue
315: !
316: ! Assemble matrix
317: !
318: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
319: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
321: !
322: ! Set operators. Here the matrix that defines the linear system
323: ! also serves as the preconditioning matrix. Since all the matrices
324: ! will have the same nonzero pattern here, we indicate this so the
325: ! linear solvers can take advantage of this.
326: !
327: call KSPSetOperators(ksp,A,A,ierr)
329: !
330: ! Set linear solver defaults for this problem (optional).
331: ! - Here we set it to use direct LU factorization for the solution
332: !
333: call KSPGetPC(ksp,pc,ierr)
334: call PCSetType(pc,PCLU,ierr)
336: !
337: ! Set runtime options, e.g.,
338: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
339: ! These options will override those specified above as long as
340: ! KSPSetFromOptions() is called _after_ any other customization
341: ! routines.
342: !
343: ! Run the program with the option -help to see all the possible
344: ! linear solver options.
345: !
346: call KSPSetFromOptions(ksp,ierr)
348: !
349: ! This allows the PETSc linear solvers to compute the solution
350: ! directly in the user's array rather than in the PETSc vector.
351: !
352: ! This is essentially a hack and not highly recommend unless you
353: ! are quite comfortable with using PETSc. In general, users should
354: ! write their entire application using PETSc vectors rather than
355: ! arrays.
356: !
357: call VecPlaceArray(x,userx,ierr)
358: call VecPlaceArray(b,userb,ierr)
360: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
361: ! Solve the linear system
362: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
364: call KSPSolve(ksp,b,x,ierr)
366: call VecResetArray(x,ierr)
367: call VecResetArray(b,ierr)
369: return
370: end
372: ! ------------------------------------------------------------------------
374: subroutine UserFinalizeLinearSolver(userctx,ierr)
376: implicit none
378: #include <petsc/finclude/petscsys.h>
379: #include <petsc/finclude/petscvec.h>
380: #include <petsc/finclude/petscmat.h>
381: #include <petsc/finclude/petscksp.h>
382: #include <petsc/finclude/petscpc.h>
384: PetscErrorCode ierr
385: PetscFortranAddr userctx(*)
387: ! Local variable declararions
389: Vec x,b
390: Mat A
391: KSP ksp
392: !
393: ! We are all done and don't need to solve any more linear systems, so
394: ! we free the work space. All PETSc objects should be destroyed when
395: ! they are no longer needed.
396: !
397: x = userctx(1)
398: b = userctx(2)
399: A = userctx(3)
400: ksp = userctx(4)
402: call VecDestroy(x,ierr)
403: call VecDestroy(b,ierr)
404: call MatDestroy(A,ierr)
405: call KSPDestroy(ksp,ierr)
407: return
408: end