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

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

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

 72: !  Input/output variables:
 73:     type(tSNES) snesIn
 74:     type(tVec) X, F
 75:     PetscErrorCode ierr
 76:     type(userctx) user

 78: !  Declarations for use with local arrays:
 79:     PetscScalar, pointer :: lx_v(:), lf_v(:)
 80:     type(tVec) 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(DMGetLocalVector(user%da, localX, ierr))
 87:     PetscCall(DMGlobalToLocalBegin(user%da, X, INSERT_VALUES, localX, ierr))
 88:     PetscCall(DMGlobalToLocalEnd(user%da, X, INSERT_VALUES, localX, ierr))

 90: !  Get a pointer to vector data.
 91: !    - VecGetArray90() returns a pointer to the data array.
 92: !    - You MUST call VecRestoreArray() when you no longer need access to
 93: !      the array.

 95:     PetscCall(VecGetArray(localX, lx_v, ierr))
 96:     PetscCall(VecGetArray(F, lf_v, ierr))

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

101: !  Restore vectors
102:     PetscCall(VecRestoreArray(localX, lx_v, ierr))
103:     PetscCall(VecRestoreArray(F, lf_v, ierr))

105: !  Insert values into global vector

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

110: !      PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr))
111: !      PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr))
112:   end subroutine formfunction
113: end module ex5f90tmodule

115: module f90moduleinterfacest
116:   use ex5f90tmodule

118:   Interface SNESSetApplicationContext
119:     Subroutine SNESSetApplicationContext(snesIn, ctx, ierr)
120: #include <petsc/finclude/petscsnes.h>
121:       use petscsnes
122:       use ex5f90tmodule
123:       type(tSNES) snesIn
124:       type(userctx) ctx
125:       PetscErrorCode ierr
126:     End Subroutine
127:   End Interface SNESSetApplicationContext

129:   Interface SNESGetApplicationContext
130:     Subroutine SNESGetApplicationContext(snesIn, ctx, ierr)
131: #include <petsc/finclude/petscsnes.h>
132:       use petscsnes
133:       use ex5f90tmodule
134:       type(tSNES) snesIn
135:       type(userctx), pointer :: ctx
136:       PetscErrorCode ierr
137:     End Subroutine
138:   End Interface SNESGetApplicationContext
139: end module f90moduleinterfacest

141: program main
142: #include <petsc/finclude/petscdmda.h>
143: #include <petsc/finclude/petscsnes.h>
144:   use petscdmda
145:   use petscsnes
146:   use ex5f90tmodule
147:   use f90moduleinterfacest
148:   implicit none
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: !                   Variable declarations
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: !
153: !  Variables:
154: !     mysnes      - nonlinear solver
155: !     x, r        - solution, residual vectors
156: !     J           - Jacobian matrix
157: !     its         - iterations for convergence
158: !     Nx, Ny      - number of preocessors in x- and y- directions
159: !     matrix_free - flag - 1 indicates matrix-free version
160: !
161:   type(tSNES) mysnes
162:   type(tVec) x, r
163:   type(tMat) J
164:   PetscErrorCode ierr
165:   PetscInt its
166:   PetscBool flg, matrix_free, set
167:   PetscInt ione, nfour
168:   PetscReal lambda_max, lambda_min
169:   type(userctx) user
170:   type(userctx), pointer:: puser
171:   type(tPetscOptions) :: options

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

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

183: !  Initialize problem parameters
184:   options%v = 0
185:   lambda_max = 6.81
186:   lambda_min = 0.0
187:   user%lambda = 6.0
188:   ione = 1
189:   nfour = 4
190:   PetscCallA(PetscOptionsGetReal(options, PETSC_NULL_CHARACTER, '-par', user%lambda, flg, ierr))
191:   PetscCheckA(user%lambda < lambda_max .and. user%lambda > lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda provided with -par is out of range')

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

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

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

204: ! This really needs only the star-type stencil, but we use the box
205: ! stencil temporarily.
206:   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_ARRAY, PETSC_NULL_INTEGER_ARRAY, user%da, ierr))
207:   PetscCallA(DMSetFromOptions(user%da, ierr))
208:   PetscCallA(DMSetUp(user%da, ierr))
209:   PetscCallA(DMDAGetInfo(user%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_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))

211: !
212: !   Visualize the distribution of the array across the processors
213: !
214: !     PetscCallA(DMView(user%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(user%da, x, ierr))
219:   PetscCallA(VecDuplicate(x, r, ierr))

221: !  Get local grid boundaries (for 2-dimensional DMDA)
222:   PetscCallA(DMDAGetCorners(user%da, user%xs, user%ys, PETSC_NULL_INTEGER, user%xm, user%ym, PETSC_NULL_INTEGER, ierr))
223:   PetscCallA(DMDAGetGhostCorners(user%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(mysnes, user, ierr))

239: !  Set function evaluation routine and vector
240:   PetscCallA(SNESSetFunction(mysnes, 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 matrix used to construct the preconditioner 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(options, PETSC_NULL_CHARACTER, '-snes_mf', matrix_free, ierr))
267:   if (.not. matrix_free) then
268:     PetscCallA(DMSetMatType(user%da, MATAIJ, ierr))
269:     PetscCallA(DMCreateMatrix(user%da, J, ierr))
270:     PetscCallA(SNESSetJacobian(mysnes, J, J, FormJacobian, user, ierr))
271:   end if

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(SNESSetFromOptions(mysnes, ierr))

279: !     Test Fortran90 wrapper for SNESSet/Get ApplicationContext()
280:   PetscCallA(PetscOptionsGetBool(options, PETSC_NULL_CHARACTER, '-test_appctx', flg, set, ierr))
281:   if (flg) then
282:     PetscCallA(SNESGetApplicationContext(mysnes, puser, ierr))
283:   end if

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

293:   PetscCallA(FormInitialGuess(mysnes, x, ierr))
294:   PetscCallA(SNESSolve(mysnes, PETSC_NULL_VEC, x, ierr))
295:   PetscCallA(SNESGetIterationNumber(mysnes, its, ierr))
296:   if (user%rank == 0) then
297:     write (6, 100) its
298:   end if
299: 100 format('Number of SNES iterations = ', i5)

301: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
302: !  Free work space.  All PETSc objects should be destroyed when they
303: !  are no longer needed.
304: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
305:   if (.not. matrix_free) PetscCallA(MatDestroy(J, ierr))
306:   PetscCallA(VecDestroy(x, ierr))
307:   PetscCallA(VecDestroy(r, ierr))
308:   PetscCallA(SNESDestroy(mysnes, ierr))
309:   PetscCallA(DMDestroy(user%da, ierr))

311:   PetscCallA(PetscFinalize(ierr))
312: end

314: ! ---------------------------------------------------------------------
315: !
316: !  FormInitialGuess - Forms initial approximation.
317: !
318: !  Input Parameters:
319: !  X - vector
320: !
321: !  Output Parameter:
322: !  X - vector
323: !
324: !  Notes:
325: !  This routine serves as a wrapper for the lower-level routine
326: !  "InitialGuessLocal", where the actual computations are
327: !  done using the standard Fortran style of treating the local
328: !  vector data as a multidimensional array over the local mesh.
329: !  This routine merely handles ghost point scatters and accesses
330: !  the local vector data via VecGetArray() and VecRestoreArray().
331: !
332: subroutine FormInitialGuess(mysnes, X, ierr)
333: #include <petsc/finclude/petscsnes.h>
334:   use petscsnes
335:   use ex5f90tmodule
336:   use f90moduleinterfacest
337: !  Input/output variables:
338:   type(tSNES) mysnes
339:   type(userctx), pointer:: puser
340:   type(tVec) X
341:   PetscErrorCode ierr

343: !  Declarations for use with local arrays:
344:   PetscScalar, pointer :: lx_v(:)

346:   ierr = 0
347:   PetscCallA(SNESGetApplicationContext(mysnes, puser, ierr))
348: !  Get a pointer to vector data.
349: !    - VecGetArray90() returns a pointer to the data array.
350: !    - You MUST call VecRestoreArray() when you no longer need access to
351: !      the array.

353:   PetscCallA(VecGetArray(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(VecRestoreArray(X, lx_v, ierr))

361: !  Insert values into global vector

363: end

365: ! ---------------------------------------------------------------------
366: !
367: !  InitialGuessLocal - Computes initial approximation, called by
368: !  the higher level routine FormInitialGuess().
369: !
370: !  Input Parameter:
371: !  x - local vector data
372: !
373: !  Output Parameters:
374: !  x - local vector data
375: !  ierr - error code
376: !
377: !  Notes:
378: !  This routine uses standard Fortran-style computations over a 2-dim array.
379: !
380: subroutine InitialGuessLocal(user, x, ierr)
381: #include <petsc/finclude/petscsys.h>
382:   use petscsys
383:   use ex5f90tmodule
384: !  Input/output variables:
385:   type(userctx) user
386:   PetscScalar x(user%xs:user%xe, user%ys:user%ye)
387:   PetscErrorCode ierr

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

394: !  Set parameters

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

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

413:     end

415: ! ---------------------------------------------------------------------
416: !
417: !  FormFunctionLocal - Computes nonlinear function, called by
418: !  the higher level routine FormFunction().
419: !
420: !  Input Parameter:
421: !  x - local vector data
422: !
423: !  Output Parameters:
424: !  f - local vector data, f(x)
425: !  ierr - error code
426: !
427: !  Notes:
428: !  This routine uses standard Fortran-style computations over a 2-dim array.
429: !
430:     subroutine FormFunctionLocal(x, f, user, ierr)
431: #include <petsc/finclude/petscsys.h>
432:       use petscsys
433:       use ex5f90tmodule
434: !  Input/output variables:
435:       type(userctx) user
436:       PetscScalar x(user%gxs:user%gxe, user%gys:user%gye)
437:       PetscScalar f(user%xs:user%xe, user%ys:user%ye)
438:       PetscErrorCode ierr

440: !  Local variables:
441:       PetscScalar two, one, hx, hy, hxdhy, hydhx, sc
442:       PetscScalar u, uxx, uyy
443:       PetscInt i, j

445:       one = 1.0
446:       two = 2.0
447:       hx = one/PetscIntToReal(user%mx - 1)
448:       hy = one/PetscIntToReal(user%my - 1)
449:       sc = hx*hy*user%lambda
450:       hxdhy = hx/hy
451:       hydhx = hy/hx

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

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

470: ! ---------------------------------------------------------------------
471: !
472: !  FormJacobian - Evaluates Jacobian matrix.
473: !
474: !  Input Parameters:
475: !  snes     - the SNES context
476: !  x        - input vector
477: !  dummy    - optional user-defined context, as set by SNESSetJacobian()
478: !             (not used here)
479: !
480: !  Output Parameters:
481: !  jac      - Jacobian matrix
482: !  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
483: !
484: !  Notes:
485: !  This routine serves as a wrapper for the lower-level routine
486: !  "FormJacobianLocal", where the actual computations are
487: !  done using the standard Fortran style of treating the local
488: !  vector data as a multidimensional array over the local mesh.
489: !  This routine merely accesses the local vector data via
490: !  VecGetArray() and VecRestoreArray().
491: !
492: !  Notes:
493: !  Due to grid point reordering with DMDAs, we must always work
494: !  with the local grid points, and then transform them to the new
495: !  global numbering with the "ltog" mapping
496: !  We cannot work directly with the global numbers for the original
497: !  uniprocessor grid!
498: !
499: !  Two methods are available for imposing this transformation
500: !  when setting matrix entries:
501: !    (A) MatSetValuesLocal(), using the local ordering (including
502: !        ghost points!)
503: !        - Set matrix entries using the local ordering
504: !          by calling MatSetValuesLocal()
505: !    (B) MatSetValues(), using the global ordering
506: !        - Use DMGetLocalToGlobalMapping() then
507: !          ISLocalToGlobalMappingGetIndices() to extract the local-to-global map
508: !        - Then apply this map explicitly yourself
509: !        - Set matrix entries using the global ordering by calling
510: !          MatSetValues()
511: !  Option (A) seems cleaner/easier in many cases, and is the procedure
512: !  used in this example.
513: !
514:         subroutine FormJacobian(mysnes, X, jac, jac_prec, user, ierr)
515: #include <petsc/finclude/petscsnes.h>
516:           use petscsnes
517:           use ex5f90tmodule
518: !  Input/output variables:
519:           type(tSNES) mysnes
520:           type(tVec) X
521:           type(tMat) jac, jac_prec
522:           type(userctx) user
523:           PetscErrorCode ierr

525: !  Declarations for use with local arrays:
526:           PetscScalar, pointer :: lx_v(:)
527:           type(tVec) localX

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

534:           PetscCallA(DMGetLocalVector(user%da, localX, ierr))
535:           PetscCallA(DMGlobalToLocalBegin(user%da, X, INSERT_VALUES, localX, ierr))
536:           PetscCallA(DMGlobalToLocalEnd(user%da, X, INSERT_VALUES, localX, ierr))

538: !  Get a pointer to vector data
539:           PetscCallA(VecGetArray(localX, lx_v, ierr))

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

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

549:           PetscCallA(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
550: !      if (jac .ne. jac_prec) then
551:           PetscCallA(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
552: !      endif
553:           PetscCallA(VecRestoreArray(localX, lx_v, ierr))
554:           PetscCallA(DMRestoreLocalVector(user%da, localX, ierr))
555:           PetscCallA(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
556: !      if (jac .ne. jac_prec) then
557:           PetscCallA(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
558: !      endif

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

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

565:         end

567: ! ---------------------------------------------------------------------
568: !
569: !  FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner,
570: !  called by the higher level routine FormJacobian().
571: !
572: !  Input Parameters:
573: !  x        - local vector data
574: !
575: !  Output Parameters:
576: !  jac_prec - Jacobian matrix used to compute the preconditioner
577: !  ierr     - error code
578: !
579: !  Notes:
580: !  This routine uses standard Fortran-style computations over a 2-dim array.
581: !
582: !  Notes:
583: !  Due to grid point reordering with DMDAs, we must always work
584: !  with the local grid points, and then transform them to the new
585: !  global numbering with the "ltog" mapping
586: !  We cannot work directly with the global numbers for the original
587: !  uniprocessor grid!
588: !
589: !  Two methods are available for imposing this transformation
590: !  when setting matrix entries:
591: !    (A) MatSetValuesLocal(), using the local ordering (including
592: !        ghost points!)
593: !        - Set matrix entries using the local ordering
594: !          by calling MatSetValuesLocal()
595: !    (B) MatSetValues(), using the global ordering
596: !        - Set matrix entries using the global ordering by calling
597: !          MatSetValues()
598: !  Option (A) seems cleaner/easier in many cases, and is the procedure
599: !  used in this example.
600: !
601:         subroutine FormJacobianLocal(x, jac_prec, user, ierr)
602: #include <petsc/finclude/petscmat.h>
603:           use petscmat
604:           use ex5f90tmodule
605: !  Input/output variables:
606:           type(userctx) user
607:           PetscScalar x(user%gxs:user%gxe, user%gys:user%gye)
608:           type(tMat) jac_prec
609:           PetscErrorCode ierr

611: !  Local variables:
612:           PetscInt row, col(5), i, j
613:           PetscInt ione, ifive
614:           PetscScalar two, one, hx, hy, hxdhy
615:           PetscScalar hydhx, sc, v(5)

617: !  Set parameters
618:           ione = 1
619:           ifive = 5
620:           one = 1.0
621:           two = 2.0
622:           hx = one/PetscIntToReal(user%mx - 1)
623:           hy = one/PetscIntToReal(user%my - 1)
624:           sc = hx*hy
625:           hxdhy = hx/hy
626:           hydhx = hy/hx

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

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

667: !/*TEST
668: !
669: !   test:
670: !      nsize: 4
671: !      args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
672: !
673: !TEST*/