Actual source code: ex13f90.F

petsc-master 2016-06-24
Report Typos and Errors
  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/(m+1)
120:       hy = 1.0/(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