Actual source code: ex1f.F90

  1: !
  2: !  Description: This example solves a nonlinear system on 1 processor with SNES.
  3: !  We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
  4: !  domain.  The command line options include:
  5: !    -par <parameter>, where <parameter> indicates the nonlinearity of the problem
  6: !       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)
  7: !    -mx <xg>, where <xg> = number of grid points in the x-direction
  8: !    -my <yg>, where <yg> = number of grid points in the y-direction
  9: !

 11: !
 12: !  --------------------------------------------------------------------------
 13: !
 14: !  Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 15: !  the partial differential equation
 16: !
 17: !          -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 18: !
 19: !  with boundary conditions
 20: !
 21: !           u = 0  for  x = 0, x = 1, y = 0, y = 1.
 22: !
 23: !  A finite difference approximation with the usual 5-point stencil
 24: !  is used to discretize the boundary value problem to obtain a nonlinear
 25: !  system of equations.
 26: !
 27: !  The parallel version of this code is snes/tutorials/ex5f.F
 28: !
 29: !  --------------------------------------------------------------------------
 30: #include <petsc/finclude/petscsnes.h>
 31: #include <petsc/finclude/petscdraw.h>
 32: module ex1fmodule
 33:   use petscsnes
 34:   implicit none
 35: contains
 36:   subroutine postcheck(snes, x, y, w, changed_y, changed_w, ctx, ierr)
 37:     SNES snes
 38:     PetscReal norm
 39:     Vec tmp, x, y, w
 40:     PetscBool changed_w, changed_y
 41:     PetscErrorCode ierr
 42:     PetscInt ctx
 43:     PetscScalar mone
 44:     MPIU_Comm comm

 46:     character(len=PETSC_MAX_PATH_LEN) :: outputString

 48:     PetscCallA(PetscObjectGetComm(snes, comm, ierr))
 49:     PetscCallA(VecDuplicate(x, tmp, ierr))
 50:     mone = -1.0
 51:     PetscCallA(VecWAXPY(tmp, mone, x, w, ierr))
 52:     PetscCallA(VecNorm(tmp, NORM_2, norm, ierr))
 53:     PetscCallA(VecDestroy(tmp, ierr))
 54:     write (outputString, *) norm
 55:     PetscCallA(PetscPrintf(comm, 'Norm of search step '//trim(outputString)//'\n', ierr))
 56:   end

 58: ! ---------------------------------------------------------------------
 59: !
 60: !  FormInitialGuess - Forms initial approximation.
 61: !
 62: !  Input Parameter:
 63: !  X - vector
 64: !
 65: !  Output Parameters:
 66: !  X - vector
 67: !  ierr - error code
 68: !
 69: !  Notes:
 70: !  This routine serves as a wrapper for the lower-level routine
 71: !  "ApplicationInitialGuess", where the actual computations are
 72: !  done using the standard Fortran style of treating the local
 73: !  vector data as a multidimensional array over the local mesh.
 74: !  This routine merely accesses the local vector data via
 75: !  VecGetArray() and VecRestoreArray().
 76: !
 77:   subroutine FormInitialGuess(X, ierr)

 79: !  Input/output variables:
 80:     Vec X
 81:     PetscErrorCode ierr

 83: !     Declarations for use with local arrays:
 84:     PetscScalar, pointer :: lx_v(:)

 86:     ierr = 0

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

 93:     PetscCallA(VecGetArray(X, lx_v, ierr))

 95: !  Compute initial guess

 97:     PetscCallA(ApplicationInitialGuess(lx_v, ierr))

 99: !  Restore vector

101:     PetscCallA(VecRestoreArray(X, lx_v, ierr))

103:   end

105: !  ApplicationInitialGuess - Computes initial approximation, called by
106: !  the higher level routine FormInitialGuess().
107: !
108: !  Input Parameter:
109: !  x - local vector data
110: !
111: !  Output Parameters:
112: !  f - local vector data, f(x)
113: !  ierr - error code
114: !
115: !  Notes:
116: !  This routine uses standard Fortran-style computations over a 2-dim array.
117: !
118:   subroutine ApplicationInitialGuess(x, ierr)

120: !  Common blocks:
121:     PetscReal lambda
122:     PetscInt mx, my
123:     PetscBool fd_coloring
124:     common/params/lambda, mx, my, fd_coloring

126: !  Input/output variables:
127:     PetscScalar x(mx, my)
128:     PetscErrorCode ierr

130: !  Local variables:
131:     PetscInt i, j
132:     PetscReal temp1, temp, hx, hy, one

134: !  Set parameters

136:     ierr = 0
137:     one = 1.0
138:     hx = one/(mx - 1)
139:     hy = one/(my - 1)
140:     temp1 = lambda/(lambda + one)

142:     do j = 1, my
143:       temp = min(j - 1, my - j)*hy
144:       do i = 1, mx
145:         if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
146:           x(i, j) = 0.0
147:         else
148:           x(i, j) = temp1*sqrt(min(min(i - 1, mx - i)*hx, temp))
149:         end if
150:       end do
151:     end do

153:   end

155: ! ---------------------------------------------------------------------
156: !
157: !  FormFunction - Evaluates nonlinear function, F(x).
158: !
159: !  Input Parameters:
160: !  snes  - the SNES context
161: !  X     - input vector
162: !  dummy - optional user-defined context, as set by SNESSetFunction()
163: !          (not used here)
164: !
165: !  Output Parameter:
166: !  F     - vector with newly computed function
167: !
168: !  Notes:
169: !  This routine serves as a wrapper for the lower-level routine
170: !  "ApplicationFunction", where the actual computations are
171: !  done using the standard Fortran style of treating the local
172: !  vector data as a multidimensional array over the local mesh.
173: !  This routine merely accesses the local vector data via
174: !  VecGetArray() and VecRestoreArray().
175: !
176:   subroutine FormFunction(snes, X, F, fdcoloring, ierr)

178: !  Input/output variables:
179:     SNES snes
180:     Vec X, F
181:     PetscErrorCode ierr
182:     MatFDColoring fdcoloring

184: !  Common blocks:
185:     PetscReal lambda
186:     PetscInt mx, my
187:     PetscBool fd_coloring
188:     common/params/lambda, mx, my, fd_coloring

190: !  Declarations for use with local arrays:
191:     PetscScalar, pointer :: lx_v(:), lf_v(:)
192:     PetscInt, pointer :: indices(:)

194: !  Get pointers to vector data.
195: !    - VecGetArray() returns a pointer to the data array.
196: !    - You MUST call VecRestoreArray() when you no longer need access to
197: !      the array.

199:     PetscCallA(VecGetArrayRead(X, lx_v, ierr))
200:     PetscCallA(VecGetArray(F, lf_v, ierr))

202: !  Compute function

204:     PetscCallA(ApplicationFunction(lx_v, lf_v, ierr))

206: !  Restore vectors

208:     PetscCallA(VecRestoreArrayRead(X, lx_v, ierr))
209:     PetscCallA(VecRestoreArray(F, lf_v, ierr))

211:     PetscCallA(PetscLogFlops(11.0d0*mx*my, ierr))
212: !
213: !     fdcoloring is in the common block and used here ONLY to test the
214: !     calls to MatFDColoringGetPerturbedColumns() and  MatFDColoringRestorePerturbedColumns()
215: !
216:     if (fd_coloring) then
217:       PetscCallA(MatFDColoringGetPerturbedColumns(fdcoloring, PETSC_NULL_INTEGER, indices, ierr))
218:       print *, 'Indices from GetPerturbedColumns'
219:       write (*, 1000) indices
220: 1000  format(50i4)
221:       PetscCallA(MatFDColoringRestorePerturbedColumns(fdcoloring, PETSC_NULL_INTEGER, indices, ierr))
222:     end if
223:   end

225: ! ---------------------------------------------------------------------
226: !
227: !  ApplicationFunction - Computes nonlinear function, called by
228: !  the higher level routine FormFunction().
229: !
230: !  Input Parameter:
231: !  x    - local vector data
232: !
233: !  Output Parameters:
234: !  f    - local vector data, f(x)
235: !  ierr - error code
236: !
237: !  Notes:
238: !  This routine uses standard Fortran-style computations over a 2-dim array.
239: !
240:   subroutine ApplicationFunction(x, f, ierr)

242: !  Common blocks:
243:     PetscReal lambda
244:     PetscInt mx, my
245:     PetscBool fd_coloring
246:     common/params/lambda, mx, my, fd_coloring

248: !  Input/output variables:
249:     PetscScalar x(mx, my), f(mx, my)
250:     PetscErrorCode ierr

252: !  Local variables:
253:     PetscScalar two, one, hx, hy
254:     PetscScalar hxdhy, hydhx, sc
255:     PetscScalar u, uxx, uyy
256:     PetscInt i, j

258:     ierr = 0
259:     one = 1.0
260:     two = 2.0
261:     hx = one/(mx - 1)
262:     hy = one/(my - 1)
263:     sc = hx*hy*lambda
264:     hxdhy = hx/hy
265:     hydhx = hy/hx

267: !  Compute function

269:     do j = 1, my
270:       do i = 1, mx
271:         if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
272:           f(i, j) = x(i, j)
273:         else
274:           u = x(i, j)
275:           uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j))
276:           uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1))
277:           f(i, j) = uxx + uyy - sc*exp(u)
278:         end if
279:       end do
280:     end do

282:   end

284: ! ---------------------------------------------------------------------
285: !
286: !  FormJacobian - Evaluates Jacobian matrix.
287: !
288: !  Input Parameters:
289: !  snes    - the SNES context
290: !  x       - input vector
291: !  dummy   - optional user-defined context, as set by SNESSetJacobian()
292: !            (not used here)
293: !
294: !  Output Parameters:
295: !  jac      - Jacobian matrix
296: !  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
297: !
298: !  Notes:
299: !  This routine serves as a wrapper for the lower-level routine
300: !  "ApplicationJacobian", where the actual computations are
301: !  done using the standard Fortran style of treating the local
302: !  vector data as a multidimensional array over the local mesh.
303: !  This routine merely accesses the local vector data via
304: !  VecGetArray() and VecRestoreArray().
305: !
306:   subroutine FormJacobian(snes, X, jac, jac_prec, dummy, ierr)

308: !  Input/output variables:
309:     SNES snes
310:     Vec X
311:     Mat jac, jac_prec
312:     PetscErrorCode ierr
313:     integer dummy

315: !  Common blocks:
316:     PetscReal lambda
317:     PetscInt mx, my
318:     PetscBool fd_coloring
319:     common/params/lambda, mx, my, fd_coloring

321: !  Declarations for use with local array:
322:     PetscScalar, pointer :: lx_v(:)

324: !  Get a pointer to vector data

326:     PetscCallA(VecGetArrayRead(X, lx_v, ierr))

328: !  Compute Jacobian entries

330:     PetscCallA(ApplicationJacobian(lx_v, jac, jac_prec, ierr))

332: !  Restore vector

334:     PetscCallA(VecRestoreArrayRead(X, lx_v, ierr))

336: !  Assemble matrix

338:     PetscCallA(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
339:     PetscCallA(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr))

341:   end

343: ! ---------------------------------------------------------------------
344: !
345: !  ApplicationJacobian - Computes Jacobian matrix, called by
346: !  the higher level routine FormJacobian().
347: !
348: !  Input Parameters:
349: !  x        - local vector data
350: !
351: !  Output Parameters:
352: !  jac      - Jacobian matrix
353: !  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
354: !  ierr     - error code
355: !
356: !  Notes:
357: !  This routine uses standard Fortran-style computations over a 2-dim array.
358: !
359:   subroutine ApplicationJacobian(x, jac, jac_prec, ierr)
360: !  Common blocks:
361:     PetscReal lambda
362:     PetscInt mx, my
363:     PetscBool fd_coloring
364:     common/params/lambda, mx, my, fd_coloring

366: !  Input/output variables:
367:     PetscScalar x(mx, my)
368:     Mat jac, jac_prec
369:     PetscErrorCode ierr

371: !  Local variables:
372:     PetscInt i, j, row(1), col(5), i1, i5
373:     PetscScalar two, one, hx, hy
374:     PetscScalar hxdhy, hydhx, sc, v(5)

376: !  Set parameters

378:     i1 = 1
379:     i5 = 5
380:     one = 1.0
381:     two = 2.0
382:     hx = one/(mx - 1)
383:     hy = one/(my - 1)
384:     sc = hx*hy
385:     hxdhy = hx/hy
386:     hydhx = hy/hx

388: !  Compute entries of the Jacobian matrix
389: !   - Here, we set all entries for a particular row at once.
390: !   - Note that MatSetValues() uses 0-based row and column numbers
391: !     in Fortran as well as in C.

393:     do j = 1, my
394:       row(1) = (j - 1)*mx - 1
395:       do i = 1, mx
396:         row(1) = row(1) + 1
397: !           boundary points
398:         if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
399:           PetscCallA(MatSetValues(jac_prec, i1, row, i1, row, [one], INSERT_VALUES, ierr))
400: !           interior grid points
401:         else
402:           v(1) = -hxdhy
403:           v(2) = -hydhx
404:           v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i, j))
405:           v(4) = -hydhx
406:           v(5) = -hxdhy
407:           col(1) = row(1) - mx
408:           col(2) = row(1) - 1
409:           col(3) = row(1)
410:           col(4) = row(1) + 1
411:           col(5) = row(1) + mx
412:           PetscCallA(MatSetValues(jac_prec, i1, row, i5, col, v, INSERT_VALUES, ierr))
413:         end if
414:       end do
415:     end do

417:   end

419: end module ex1fmodule
420: program main
421:   use petscdraw
422:   use petscsnes
423:   use ex1fmodule
424:   implicit none
425: !
426: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
427: !                   Variable declarations
428: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
429: !
430: !  Variables:
431: !     snes        - nonlinear solver
432: !     x, r        - solution, residual vectors
433: !     J           - Jacobian matrix
434: !     its         - iterations for convergence
435: !     matrix_free - flag - 1 indicates matrix-free version
436: !     lambda      - nonlinearity parameter
437: !     draw        - drawing context
438: !
439:   SNES snes
440:   MatColoring mc
441:   Vec x, r
442:   PetscDraw draw
443:   Mat J
444:   PetscBool matrix_free, flg, fd_coloring
445:   PetscErrorCode ierr
446:   PetscInt its, N, mx, my, i5
447:   PetscMPIInt size, rank
448:   PetscReal lambda_max, lambda_min, lambda
449:   MatFDColoring fdcoloring
450:   ISColoring iscoloring
451:   PetscBool pc
452:   integer4 imx, imy
453:   character(len=PETSC_MAX_PATH_LEN) :: outputString
454:   PetscScalar, pointer :: lx_v(:)
455:   integer4 xl, yl, width, height

457: !  Store parameters in common block

459:   common/params/lambda, mx, my, fd_coloring

461: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
462: !  Initialize program
463: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

465:   PetscCallA(PetscInitialize(ierr))
466:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
467:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))

469:   PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only')

471: !  Initialize problem parameters
472:   i5 = 5
473:   lambda_max = 6.81
474:   lambda_min = 0.0
475:   lambda = 6.0
476:   mx = 4
477:   my = 4
478:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr))
479:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr))
480:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', lambda, flg, ierr))
481:   PetscCheckA(lambda < lambda_max .and. lambda > lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda out of range ')
482:   N = mx*my
483:   pc = PETSC_FALSE
484:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-pc', pc, PETSC_NULL_BOOL, ierr))

486: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
487: !  Create nonlinear solver context
488: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

490:   PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))

492:   if (pc .eqv. PETSC_TRUE) then
493:     PetscCallA(SNESSetType(snes, SNESNEWTONTR, ierr))
494:     PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck, snes, ierr))
495:   end if

497: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
498: !  Create vector data structures; set function evaluation routine
499: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

501:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
502:   PetscCallA(VecSetSizes(x, PETSC_DECIDE, N, ierr))
503:   PetscCallA(VecSetFromOptions(x, ierr))
504:   PetscCallA(VecDuplicate(x, r, ierr))

506: !  Set function evaluation routine and vector.  Whenever the nonlinear
507: !  solver needs to evaluate the nonlinear function, it will call this
508: !  routine.
509: !   - Note that the final routine argument is the user-defined
510: !     context that provides application-specific data for the
511: !     function evaluation routine.

513:   PetscCallA(SNESSetFunction(snes, r, FormFunction, fdcoloring, ierr))

515: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
516: !  Create matrix data structure; set Jacobian evaluation routine
517: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

519: !  Create matrix. Here we only approximately preallocate storage space
520: !  for the Jacobian.  See the users manual for a discussion of better
521: !  techniques for preallocating matrix memory.

523:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-snes_mf', matrix_free, ierr))
524:   if (.not. matrix_free) then
525:     PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD, N, N, i5, PETSC_NULL_INTEGER_ARRAY, J, ierr))
526:   end if

528: !
529: !     This option will cause the Jacobian to be computed via finite differences
530: !    efficiently using a coloring of the columns of the matrix.
531: !
532:   fd_coloring = .false.
533:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-snes_fd_coloring', fd_coloring, ierr))
534:   if (fd_coloring) then

536: !
537: !      This initializes the nonzero structure of the Jacobian. This is artificial
538: !      because clearly if we had a routine to compute the Jacobian we won't need
539: !      to use finite differences.
540: !
541:     PetscCallA(FormJacobian(snes, x, J, J, 0, ierr))
542: !
543: !       Color the matrix, i.e. determine groups of columns that share no common
544: !      rows. These columns in the Jacobian can all be computed simultaneously.
545: !
546:     PetscCallA(MatColoringCreate(J, mc, ierr))
547:     PetscCallA(MatColoringSetType(mc, MATCOLORINGNATURAL, ierr))
548:     PetscCallA(MatColoringSetFromOptions(mc, ierr))
549:     PetscCallA(MatColoringApply(mc, iscoloring, ierr))
550:     PetscCallA(MatColoringDestroy(mc, ierr))
551: !
552: !       Create the data structure that SNESComputeJacobianDefaultColor() uses
553: !       to compute the actual Jacobians via finite differences.
554: !
555:     PetscCallA(MatFDColoringCreate(J, iscoloring, fdcoloring, ierr))
556:     PetscCallA(MatFDColoringSetFunction(fdcoloring, FormFunction, fdcoloring, ierr))
557:     PetscCallA(MatFDColoringSetFromOptions(fdcoloring, ierr))
558:     PetscCallA(MatFDColoringSetUp(J, iscoloring, fdcoloring, ierr))
559: !
560: !        Tell SNES to use the routine SNESComputeJacobianDefaultColor()
561: !      to compute Jacobians.
562: !
563:     PetscCallA(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, fdcoloring, ierr))
564:     PetscCallA(ISColoringDestroy(iscoloring, ierr))

566:   else if (.not. matrix_free) then

568: !  Set Jacobian matrix data structure and default Jacobian evaluation
569: !  routine.  Whenever the nonlinear solver needs to compute the
570: !  Jacobian matrix, it will call this routine.
571: !   - Note that the final routine argument is the user-defined
572: !     context that provides application-specific data for the
573: !     Jacobian evaluation routine.
574: !   - The user can override with:
575: !      -snes_fd : default finite differencing approximation of Jacobian
576: !      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
577: !                 (unless user explicitly sets preconditioner)
578: !      -snes_mf_operator : form matrix from which to construct the preconditioner as set by the user,
579: !                          but use matrix-free approx for Jacobian-vector
580: !                          products within Newton-Krylov method
581: !
582:     PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr))
583:   end if

585: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
586: !  Customize nonlinear solver; set runtime options
587: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

589: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)

591:   PetscCallA(SNESSetFromOptions(snes, ierr))

593: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
594: !  Evaluate initial guess; then solve nonlinear system.
595: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

597: !  Note: The user should initialize the vector, x, with the initial guess
598: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
599: !  to employ an initial guess of zero, the user should explicitly set
600: !  this vector to zero by calling VecSet().

602:   PetscCallA(FormInitialGuess(x, ierr))
603:   PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
604:   PetscCallA(SNESGetIterationNumber(snes, its, ierr))
605:   write (outputString, *) its
606:   PetscCallA(PetscPrintf(PETSC_COMM_WORLD, 'Number of SNES iterations = '//trim(outputString)//'\n', ierr))

608: !  PetscDraw contour plot of solution

610:   xl = 300
611:   yl = 0
612:   width = 300
613:   height = 300
614:   PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, 'Solution', xl, yl, width, height, draw, ierr))
615:   PetscCallA(PetscDrawSetFromOptions(draw, ierr))

617:   PetscCallA(VecGetArrayRead(x, lx_v, ierr))
618:   imx = int(mx, kind=kind(imx))
619:   imy = int(my, kind=kind(imy))
620:   PetscCallA(PetscDrawTensorContour(draw, imx, imy, PETSC_NULL_REAL_ARRAY, PETSC_NULL_REAL_ARRAY, lx_v, ierr))
621:   PetscCallA(VecRestoreArrayRead(x, lx_v, ierr))

623: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
624: !  Free work space.  All PETSc objects should be destroyed when they
625: !  are no longer needed.
626: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

628:   if (.not. matrix_free) PetscCallA(MatDestroy(J, ierr))
629:   if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring, ierr))

631:   PetscCallA(VecDestroy(x, ierr))
632:   PetscCallA(VecDestroy(r, ierr))
633:   PetscCallA(SNESDestroy(snes, ierr))
634:   PetscCallA(PetscDrawDestroy(draw, ierr))
635:   PetscCallA(PetscFinalize(ierr))
636: end
637: !
638: !/*TEST
639: !
640: !   build:
641: !      requires: !single !complex
642: !
643: !   test:
644: !      args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
645: !
646: !   test:
647: !      suffix: 2
648: !      args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
649: !
650: !   test:
651: !      suffix: 3
652: !      args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
653: !      filter: sort -b
654: !      filter_output: sort -b
655: !
656: !   test:
657: !     suffix: 4
658: !     args: -pc -par 6.807 -nox
659: !
660: !TEST*/