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