Actual source code: ex5f90.F90

  1: !
  2: !  Description: Solves a nonlinear system in parallel with SNES.
  3: !  We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
  4: !  domain, using distributed arrays (DMDAs) to partition the parallel grid.
  5: !  The command line options include:
  6: !    -par <parameter>, where <parameter> indicates the nonlinearity of the problem
  7: !       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)
  8: !

 10: !
 11: !  --------------------------------------------------------------------------
 12: !
 13: !  Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 14: !  the partial differential equation
 15: !
 16: !          -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 17: !
 18: !  with boundary conditions
 19: !
 20: !           u = 0  for  x = 0, x = 1, y = 0, y = 1.
 21: !
 22: !  A finite difference approximation with the usual 5-point stencil
 23: !  is used to discretize the boundary value problem to obtain a nonlinear
 24: !  system of equations.
 25: !
 26: !  The uniprocessor version of this code is snes/tutorials/ex4f.F
 27: !
 28: !  --------------------------------------------------------------------------
 29: !  The following define must be used before including any PETSc include files
 30: !  into a module or interface. This is because they can't handle declarations
 31: !  in them
 32: !

 34:       module ex5f90module
 35:       use petscsnes
 36:       use petscdmda
 37: #include <petsc/finclude/petscsnes.h>
 38:       type userctx
 39:         PetscInt xs,xe,xm,gxs,gxe,gxm
 40:         PetscInt ys,ye,ym,gys,gye,gym
 41:         PetscInt mx,my
 42:         PetscMPIInt rank
 43:         PetscReal lambda
 44:       end type userctx

 46:       contains
 47: ! ---------------------------------------------------------------------
 48: !
 49: !  FormFunction - Evaluates nonlinear function, F(x).
 50: !
 51: !  Input Parameters:
 52: !  snes - the SNES context
 53: !  X - input vector
 54: !  dummy - optional user-defined context, as set by SNESSetFunction()
 55: !          (not used here)
 56: !
 57: !  Output Parameter:
 58: !  F - function vector
 59: !
 60: !  Notes:
 61: !  This routine serves as a wrapper for the lower-level routine
 62: !  "FormFunctionLocal", where the actual computations are
 63: !  done using the standard Fortran style of treating the local
 64: !  vector data as a multidimensional array over the local mesh.
 65: !  This routine merely handles ghost point scatters and accesses
 66: !  the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
 67: !
 68:       subroutine FormFunction(snes,X,F,user,ierr)
 69:       implicit none

 71: !  Input/output variables:
 72:       SNES           snes
 73:       Vec            X,F
 74:       PetscErrorCode ierr
 75:       type (userctx) user
 76:       DM             da

 78: !  Declarations for use with local arrays:
 79:       PetscScalar,pointer :: lx_v(:),lf_v(:)
 80:       Vec            localX

 82: !  Scatter ghost points to local vector, using the 2-step process
 83: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
 84: !  By placing code between these two statements, computations can
 85: !  be done while messages are in transition.
 86:       PetscCall(SNESGetDM(snes,da,ierr))
 87:       PetscCall(DMGetLocalVector(da,localX,ierr))
 88:       PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr))
 89:       PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr))

 91: !  Get a pointer to vector data.
 92: !    - For default PETSc vectors, VecGetArray90() returns a pointer to
 93: !      the data array. Otherwise, the routine is implementation dependent.
 94: !    - You MUST call VecRestoreArrayF90() when you no longer need access to
 95: !      the array.
 96: !    - Note that the interface to VecGetArrayF90() differs from VecGetArray(),
 97: !      and is useable from Fortran-90 Only.

 99:       PetscCall(VecGetArrayF90(localX,lx_v,ierr))
100:       PetscCall(VecGetArrayF90(F,lf_v,ierr))

102: !  Compute function over the locally owned part of the grid
103:       PetscCall(FormFunctionLocal(lx_v,lf_v,user,ierr))

105: !  Restore vectors
106:       PetscCall(VecRestoreArrayF90(localX,lx_v,ierr))
107:       PetscCall(VecRestoreArrayF90(F,lf_v,ierr))

109: !  Insert values into global vector

111:       PetscCall(DMRestoreLocalVector(da,localX,ierr))
112:       PetscCall(PetscLogFlops(11.0d0*user%ym*user%xm,ierr))

114: !      PetscCallA(VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr))
115: !      PetscCallA(VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr))
116:       return
117:       end subroutine formfunction
118:       end module ex5f90module

120:       module ex5f90moduleinterfaces
121:         use ex5f90module

123:       Interface SNESSetApplicationContext
124:         Subroutine SNESSetApplicationContext(snes,ctx,ierr)
125:         use ex5f90module
126:           SNES snes
127:           type(userctx) ctx
128:           PetscErrorCode ierr
129:         End Subroutine
130:       End Interface SNESSetApplicationContext

132:       Interface SNESGetApplicationContext
133:         Subroutine SNESGetApplicationContext(snes,ctx,ierr)
134:         use ex5f90module
135:           SNES snes
136:           type(userctx), pointer :: ctx
137:           PetscErrorCode ierr
138:         End Subroutine
139:       End Interface SNESGetApplicationContext
140:       end module ex5f90moduleinterfaces

142:       program main
143:       use ex5f90module
144:       use ex5f90moduleinterfaces
145:       implicit none
146: !

148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: !                   Variable declarations
150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: !
152: !  Variables:
153: !     snes        - nonlinear solver
154: !     x, r        - solution, residual vectors
155: !     J           - Jacobian matrix
156: !     its         - iterations for convergence
157: !     Nx, Ny      - number of preocessors in x- and y- directions
158: !     matrix_free - flag - 1 indicates matrix-free version
159: !
160:       SNES             snes
161:       Vec              x,r
162:       Mat              J
163:       PetscErrorCode   ierr
164:       PetscInt         its
165:       PetscBool        flg,matrix_free
166:       PetscInt         ione,nfour
167:       PetscReal lambda_max,lambda_min
168:       type (userctx)   user
169:       DM               da

171: !  Note: Any user-defined Fortran routines (such as FormJacobian)
172: !  MUST be declared as external.
173:       external FormInitialGuess,FormJacobian

175: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: !  Initialize program
177: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:       PetscCallA(PetscInitialize(ierr))
179:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,user%rank,ierr))

181: !  Initialize problem parameters
182:       lambda_max  = 6.81
183:       lambda_min  = 0.0
184:       user%lambda = 6.0
185:       ione = 1
186:       nfour = 4
187:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',user%lambda,flg,ierr))
188:       if (user%lambda .ge. lambda_max .or. user%lambda .le. lambda_min) then
189:          SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range ')
190:       endif

192: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193: !  Create nonlinear solver context
194: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195:       PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))

197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: !  Create vector data structures; set function evaluation routine
199: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

201: !  Create distributed array (DMDA) to manage parallel grid and vectors

203: ! This really needs only the star-type stencil, but we use the box
204: ! stencil temporarily.
205:       PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nfour,nfour,PETSC_DECIDE,PETSC_DECIDE,ione,ione,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr))
206:       PetscCallA(DMSetFromOptions(da,ierr))
207:       PetscCallA(DMSetUp(da,ierr))

209:       PetscCallA(DMDAGetInfo(da,PETSC_NULL_INTEGER,user%mx,user%my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))

211: !
212: !   Visualize the distribution of the array across the processors
213: !
214: !     PetscCallA(DMView(da,PETSC_VIEWER_DRAW_WORLD,ierr))

216: !  Extract global and local vectors from DMDA; then duplicate for remaining
217: !  vectors that are the same types
218:       PetscCallA(DMCreateGlobalVector(da,x,ierr))
219:       PetscCallA(VecDuplicate(x,r,ierr))

221: !  Get local grid boundaries (for 2-dimensional DMDA)
222:       PetscCallA(DMDAGetCorners(da,user%xs,user%ys,PETSC_NULL_INTEGER,user%xm,user%ym,PETSC_NULL_INTEGER,ierr))
223:       PetscCallA(DMDAGetGhostCorners(da,user%gxs,user%gys,PETSC_NULL_INTEGER,user%gxm,user%gym,PETSC_NULL_INTEGER,ierr))

225: !  Here we shift the starting indices up by one so that we can easily
226: !  use the Fortran convention of 1-based indices (rather 0-based indices).
227:       user%xs  = user%xs+1
228:       user%ys  = user%ys+1
229:       user%gxs = user%gxs+1
230:       user%gys = user%gys+1

232:       user%ye  = user%ys+user%ym-1
233:       user%xe  = user%xs+user%xm-1
234:       user%gye = user%gys+user%gym-1
235:       user%gxe = user%gxs+user%gxm-1

237:       PetscCallA(SNESSetApplicationContext(snes,user,ierr))

239: !  Set function evaluation routine and vector
240:       PetscCallA(SNESSetFunction(snes,r,FormFunction,user,ierr))

242: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: !  Create matrix data structure; set Jacobian evaluation routine
244: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

246: !  Set Jacobian matrix data structure and default Jacobian evaluation
247: !  routine. User can override with:
248: !     -snes_fd : default finite differencing approximation of Jacobian
249: !     -snes_mf : matrix-free Newton-Krylov method with no preconditioning
250: !                (unless user explicitly sets preconditioner)
251: !     -snes_mf_operator : form preconditioning matrix as set by the user,
252: !                         but use matrix-free approx for Jacobian-vector
253: !                         products within Newton-Krylov method
254: !
255: !  Note:  For the parallel case, vectors and matrices MUST be partitioned
256: !     accordingly.  When using distributed arrays (DMDAs) to create vectors,
257: !     the DMDAs determine the problem partitioning.  We must explicitly
258: !     specify the local matrix dimensions upon its creation for compatibility
259: !     with the vector distribution.  Thus, the generic MatCreate() routine
260: !     is NOT sufficient when working with distributed arrays.
261: !
262: !     Note: Here we only approximately preallocate storage space for the
263: !     Jacobian.  See the users manual for a discussion of better techniques
264: !     for preallocating matrix memory.

266:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr))
267:       if (.not. matrix_free) then
268:         PetscCallA(DMSetMatType(da,MATAIJ,ierr))
269:         PetscCallA(DMCreateMatrix(da,J,ierr))
270:         PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,user,ierr))
271:       endif

273: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274: !  Customize nonlinear solver; set runtime options
275: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
277:       PetscCallA(SNESSetDM(snes,da,ierr))
278:       PetscCallA(SNESSetFromOptions(snes,ierr))

280: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281: !  Evaluate initial guess; then solve nonlinear system.
282: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
283: !  Note: The user should initialize the vector, x, with the initial guess
284: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
285: !  to employ an initial guess of zero, the user should explicitly set
286: !  this vector to zero by calling VecSet().

288:       PetscCallA(FormInitialGuess(snes,x,ierr))
289:       PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
290:       PetscCallA(SNESGetIterationNumber(snes,its,ierr))
291:       if (user%rank .eq. 0) then
292:          write(6,100) its
293:       endif
294:   100 format('Number of SNES iterations = ',i5)

296: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
297: !  Free work space.  All PETSc objects should be destroyed when they
298: !  are no longer needed.
299: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
300:       if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr))
301:       PetscCallA(VecDestroy(x,ierr))
302:       PetscCallA(VecDestroy(r,ierr))
303:       PetscCallA(SNESDestroy(snes,ierr))
304:       PetscCallA(DMDestroy(da,ierr))

306:       PetscCallA(PetscFinalize(ierr))
307:       end

309: ! ---------------------------------------------------------------------
310: !
311: !  FormInitialGuess - Forms initial approximation.
312: !
313: !  Input Parameters:
314: !  X - vector
315: !
316: !  Output Parameter:
317: !  X - vector
318: !
319: !  Notes:
320: !  This routine serves as a wrapper for the lower-level routine
321: !  "InitialGuessLocal", where the actual computations are
322: !  done using the standard Fortran style of treating the local
323: !  vector data as a multidimensional array over the local mesh.
324: !  This routine merely handles ghost point scatters and accesses
325: !  the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
326: !
327:       subroutine FormInitialGuess(snes,X,ierr)
328:       use ex5f90module
329:       use ex5f90moduleinterfaces
330:       implicit none

332: !  Input/output variables:
333:       SNES           snes
334:       type(userctx), pointer:: puser
335:       Vec            X
336:       PetscErrorCode ierr
337:       DM             da

339: !  Declarations for use with local arrays:
340:       PetscScalar,pointer :: lx_v(:)

342:       0
343:       PetscCallA(SNESGetDM(snes,da,ierr))
344:       PetscCallA(SNESGetApplicationContext(snes,puser,ierr))
345: !  Get a pointer to vector data.
346: !    - For default PETSc vectors, VecGetArray90() returns a pointer to
347: !      the data array. Otherwise, the routine is implementation dependent.
348: !    - You MUST call VecRestoreArrayF90() when you no longer need access to
349: !      the array.
350: !    - Note that the interface to VecGetArrayF90() differs from VecGetArray(),
351: !      and is useable from Fortran-90 Only.

353:       PetscCallA(VecGetArrayF90(X,lx_v,ierr))

355: !  Compute initial guess over the locally owned part of the grid
356:       PetscCallA(InitialGuessLocal(puser,lx_v,ierr))

358: !  Restore vector
359:       PetscCallA(VecRestoreArrayF90(X,lx_v,ierr))

361: !  Insert values into global vector

363:       return
364:       end

366: ! ---------------------------------------------------------------------
367: !
368: !  InitialGuessLocal - Computes initial approximation, called by
369: !  the higher level routine FormInitialGuess().
370: !
371: !  Input Parameter:
372: !  x - local vector data
373: !
374: !  Output Parameters:
375: !  x - local vector data
376: !  ierr - error code
377: !
378: !  Notes:
379: !  This routine uses standard Fortran-style computations over a 2-dim array.
380: !
381:       subroutine InitialGuessLocal(user,x,ierr)
382:       use ex5f90module
383:       implicit none

385: !  Input/output variables:
386:       type (userctx)         user
387:       PetscScalar  x(user%xs:user%xe,user%ys:user%ye)
388:       PetscErrorCode ierr

390: !  Local variables:
391:       PetscInt  i,j
392:       PetscReal   temp1,temp,hx,hy
393:       PetscReal   one

395: !  Set parameters

397:       0
398:       one    = 1.0
399:       hx     = one/(user%mx-1)
400:       hy     = one/(user%my-1)
401:       temp1  = user%lambda/(user%lambda + one)

403:       do 20 j=user%ys,user%ye
404:          temp = min(j-1,user%my-j)*hy
405:          do 10 i=user%xs,user%xe
406:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then
407:               x(i,j) = 0.0
408:             else
409:               x(i,j) = temp1 * sqrt(min(hx*min(i-1,user%mx-i),temp))
410:             endif
411:  10      continue
412:  20   continue

414:       return
415:       end

417: ! ---------------------------------------------------------------------
418: !
419: !  FormFunctionLocal - Computes nonlinear function, called by
420: !  the higher level routine FormFunction().
421: !
422: !  Input Parameter:
423: !  x - local vector data
424: !
425: !  Output Parameters:
426: !  f - local vector data, f(x)
427: !  ierr - error code
428: !
429: !  Notes:
430: !  This routine uses standard Fortran-style computations over a 2-dim array.
431: !
432:       subroutine FormFunctionLocal(x,f,user,ierr)
433:       use ex5f90module

435:       implicit none

437: !  Input/output variables:
438:       type (userctx) user
439:       PetscScalar  x(user%gxs:user%gxe,user%gys:user%gye)
440:       PetscScalar  f(user%xs:user%xe,user%ys:user%ye)
441:       PetscErrorCode ierr

443: !  Local variables:
444:       PetscScalar two,one,hx,hy,hxdhy,hydhx,sc
445:       PetscScalar u,uxx,uyy
446:       PetscInt  i,j

448:       one    = 1.0
449:       two    = 2.0
450:       hx     = one/(user%mx-1)
451:       hy     = one/(user%my-1)
452:       sc     = hx*hy*user%lambda
453:       hxdhy  = hx/hy
454:       hydhx  = hy/hx

456: !  Compute function over the locally owned part of the grid

458:       do 20 j=user%ys,user%ye
459:          do 10 i=user%xs,user%xe
460:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then
461:                f(i,j) = x(i,j)
462:             else
463:                u = x(i,j)
464:                uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
465:                uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
466:                f(i,j) = uxx + uyy - sc*exp(u)
467:             endif
468:  10      continue
469:  20   continue

471:       return
472:       end

474: ! ---------------------------------------------------------------------
475: !
476: !  FormJacobian - Evaluates Jacobian matrix.
477: !
478: !  Input Parameters:
479: !  snes     - the SNES context
480: !  x        - input vector
481: !  dummy    - optional user-defined context, as set by SNESSetJacobian()
482: !             (not used here)
483: !
484: !  Output Parameters:
485: !  jac      - Jacobian matrix
486: !  jac_prec - optionally different preconditioning matrix (not used here)
487: !  flag     - flag indicating matrix structure
488: !
489: !  Notes:
490: !  This routine serves as a wrapper for the lower-level routine
491: !  "FormJacobianLocal", where the actual computations are
492: !  done using the standard Fortran style of treating the local
493: !  vector data as a multidimensional array over the local mesh.
494: !  This routine merely accesses the local vector data via
495: !  VecGetArrayF90() and VecRestoreArrayF90().
496: !
497: !  Notes:
498: !  Due to grid point reordering with DMDAs, we must always work
499: !  with the local grid points, and then transform them to the new
500: !  global numbering with the "ltog" mapping
501: !  We cannot work directly with the global numbers for the original
502: !  uniprocessor grid!
503: !
504: !  Two methods are available for imposing this transformation
505: !  when setting matrix entries:
506: !    (A) MatSetValuesLocal(), using the local ordering (including
507: !        ghost points!)
508: !        - Set matrix entries using the local ordering
509: !          by calling MatSetValuesLocal()
510: !    (B) MatSetValues(), using the global ordering

512: !        - Set matrix entries using the global ordering by calling
513: !          MatSetValues()
514: !  Option (A) seems cleaner/easier in many cases, and is the procedure
515: !  used in this example.
516: !
517:       subroutine FormJacobian(snes,X,jac,jac_prec,user,ierr)
518:       use ex5f90module
519:       implicit none

521: !  Input/output variables:
522:       SNES         snes
523:       Vec          X
524:       Mat          jac,jac_prec
525:       type(userctx)  user
526:       PetscErrorCode ierr
527:       DM             da

529: !  Declarations for use with local arrays:
530:       PetscScalar,pointer :: lx_v(:)
531:       Vec            localX

533: !  Scatter ghost points to local vector, using the 2-step process
534: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd()
535: !  Computations can be done while messages are in transition,
536: !  by placing code between these two statements.

538:       PetscCallA(SNESGetDM(snes,da,ierr))
539:       PetscCallA(DMGetLocalVector(da,localX,ierr))
540:       PetscCallA(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr))
541:       PetscCallA(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr))

543: !  Get a pointer to vector data
544:       PetscCallA(VecGetArrayF90(localX,lx_v,ierr))

546: !  Compute entries for the locally owned part of the Jacobian preconditioner.
547:       PetscCallA(FormJacobianLocal(lx_v,jac_prec,user,ierr))

549: !  Assemble matrix, using the 2-step process:
550: !     MatAssemblyBegin(), MatAssemblyEnd()
551: !  Computations can be done while messages are in transition,
552: !  by placing code between these two statements.

554:       PetscCallA(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
555:       if (jac .ne. jac_prec) then
556:          PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
557:       endif
558:       PetscCallA(VecRestoreArrayF90(localX,lx_v,ierr))
559:       PetscCallA(DMRestoreLocalVector(da,localX,ierr))
560:       PetscCallA(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
561:       if (jac .ne. jac_prec) then
562:         PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
563:       endif

565: !  Tell the matrix we will never add a new nonzero location to the
566: !  matrix. If we do it will generate an error.

568:       PetscCallA(MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))

570:       return
571:       end

573: ! ---------------------------------------------------------------------
574: !
575: !  FormJacobianLocal - Computes Jacobian preconditioner matrix,
576: !  called by the higher level routine FormJacobian().
577: !
578: !  Input Parameters:
579: !  x        - local vector data
580: !
581: !  Output Parameters:
582: !  jac_prec - Jacobian preconditioner matrix
583: !  ierr     - error code
584: !
585: !  Notes:
586: !  This routine uses standard Fortran-style computations over a 2-dim array.
587: !
588: !  Notes:
589: !  Due to grid point reordering with DMDAs, we must always work
590: !  with the local grid points, and then transform them to the new
591: !  global numbering with the "ltog" mapping
592: !  We cannot work directly with the global numbers for the original
593: !  uniprocessor grid!
594: !
595: !  Two methods are available for imposing this transformation
596: !  when setting matrix entries:
597: !    (A) MatSetValuesLocal(), using the local ordering (including
598: !        ghost points!)
599: !        - Set matrix entries using the local ordering
600: !          by calling MatSetValuesLocal()
601: !    (B) MatSetValues(), using the global ordering
602: !        - Then apply this map explicitly yourself
603: !        - Set matrix entries using the global ordering by calling
604: !          MatSetValues()
605: !  Option (A) seems cleaner/easier in many cases, and is the procedure
606: !  used in this example.
607: !
608:       subroutine FormJacobianLocal(x,jac_prec,user,ierr)
609:       use ex5f90module
610:       implicit none

612: !  Input/output variables:
613:       type (userctx) user
614:       PetscScalar    x(user%gxs:user%gxe,user%gys:user%gye)
615:       Mat            jac_prec
616:       PetscErrorCode ierr

618: !  Local variables:
619:       PetscInt    row,col(5),i,j
620:       PetscInt    ione,ifive
621:       PetscScalar two,one,hx,hy,hxdhy
622:       PetscScalar hydhx,sc,v(5)

624: !  Set parameters
625:       ione   = 1
626:       ifive  = 5
627:       one    = 1.0
628:       two    = 2.0
629:       hx     = one/(user%mx-1)
630:       hy     = one/(user%my-1)
631:       sc     = hx*hy
632:       hxdhy  = hx/hy
633:       hydhx  = hy/hx

635: !  Compute entries for the locally owned part of the Jacobian.
636: !   - Currently, all PETSc parallel matrix formats are partitioned by
637: !     contiguous chunks of rows across the processors.
638: !   - Each processor needs to insert only elements that it owns
639: !     locally (but any non-local elements will be sent to the
640: !     appropriate processor during matrix assembly).
641: !   - Here, we set all entries for a particular row at once.
642: !   - We can set matrix entries either using either
643: !     MatSetValuesLocal() or MatSetValues(), as discussed above.
644: !   - Note that MatSetValues() uses 0-based row and column numbers
645: !     in Fortran as well as in C.

647:       do 20 j=user%ys,user%ye
648:          row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1
649:          do 10 i=user%xs,user%xe
650:             row = row + 1
651: !           boundary points
652:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then
653:                col(1) = row
654:                v(1)   = one
655:                PetscCallA(MatSetValuesLocal(jac_prec,ione,row,ione,col,v,INSERT_VALUES,ierr))
656: !           interior grid points
657:             else
658:                v(1) = -hxdhy
659:                v(2) = -hydhx
660:                v(3) = two*(hydhx + hxdhy) - sc*user%lambda*exp(x(i,j))
661:                v(4) = -hydhx
662:                v(5) = -hxdhy
663:                col(1) = row - user%gxm
664:                col(2) = row - 1
665:                col(3) = row
666:                col(4) = row + 1
667:                col(5) = row + user%gxm
668:                PetscCallA(MatSetValuesLocal(jac_prec,ione,row,ifive,col,v,INSERT_VALUES,ierr))
669:             endif
670:  10      continue
671:  20   continue

673:       return
674:       end

676: !
677: !/*TEST
678: !
679: !   test:
680: !      nsize: 4
681: !      args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
682: !      requires: !single
683: !
684: !   test:
685: !      suffix: 2
686: !      nsize: 4
687: !      args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
688: !      requires: !single
689: !
690: !   test:
691: !      suffix: 3
692: !      nsize: 3
693: !      args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
694: !      requires: !single
695: !
696: !   test:
697: !      suffix: 4
698: !      nsize: 3
699: !      args: -snes_mf_operator -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
700: !      requires: !single
701: !
702: !   test:
703: !      suffix: 5
704: !      requires: !single
705: !
706: !TEST*/