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: !
 10: !!/*T
 11: !  Concepts: SNES^sequential Bratu example
 12: !  Processors: 1
 13: !T*/

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

 46:       call VecDuplicate(x,tmp,ierr)
 47:       mone = -1.0
 48:       call VecWAXPY(tmp,mone,x,w,ierr)
 49:       call VecNorm(tmp,NORM_2,norm,ierr)
 50:       call VecDestroy(tmp,ierr)
 51:       print*, 'Norm of search step ',norm
 52:       changed_y = PETSC_FALSE
 53:       changed_w = PETSC_FALSE
 54:       return
 55:       end

 57:       program main
 58: #include <petsc/finclude/petscdraw.h>
 59:       use petscsnes
 60:       implicit none
 61:       interface SNESSetJacobian
 62:       subroutine SNESSetJacobian1(a,b,c,d,e,z)
 63:        use petscsnes
 64:        SNES a
 65:        Mat b
 66:        Mat c
 67:        external d
 68:        MatFDColoring e
 69:        PetscErrorCode z
 70:       end subroutine
 71:       subroutine SNESSetJacobian2(a,b,c,d,e,z)
 72:        use petscsnes
 73:        SNES a
 74:        Mat b
 75:        Mat c
 76:        external d
 77:        integer e
 78:        PetscErrorCode z
 79:       end subroutine
 80:       end interface
 81: !
 82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83: !                   Variable declarations
 84: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85: !
 86: !  Variables:
 87: !     snes        - nonlinear solver
 88: !     x, r        - solution, residual vectors
 89: !     J           - Jacobian matrix
 90: !     its         - iterations for convergence
 91: !     matrix_free - flag - 1 indicates matrix-free version
 92: !     lambda      - nonlinearity parameter
 93: !     draw        - drawing context
 94: !
 95:       SNES               snes
 96:       MatColoring        mc
 97:       Vec                x,r
 98:       PetscDraw               draw
 99:       Mat                J
100:       PetscBool  matrix_free,flg,fd_coloring
101:       PetscErrorCode ierr
102:       PetscInt   its,N, mx,my,i5
103:       PetscMPIInt size,rank
104:       PetscReal   lambda_max,lambda_min,lambda
105:       MatFDColoring      fdcoloring
106:       ISColoring         iscoloring
107:       PetscBool          pc
108:       external           postcheck

110:       PetscScalar        lx_v(0:1)
111:       PetscOffset        lx_i

113: !  Store parameters in common block

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

117: !  Note: Any user-defined Fortran routines (such as FormJacobian)
118: !  MUST be declared as external.

120:       external FormFunction,FormInitialGuess,FormJacobian

122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: !  Initialize program
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

126:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
127:       if (ierr .ne. 0) then
128:         print*,'Unable to initialize PETSc'
129:         stop
130:       endif
131:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
132:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

134:       if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif

136: !  Initialize problem parameters
137:       i5 = 5
138:       lambda_max = 6.81
139:       lambda_min = 0.0
140:       lambda     = 6.0
141:       mx         = 4
142:       my         = 4
143:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
144:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
145:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)
146:       if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range '); endif
147:       N  = mx*my
148:       pc = PETSC_FALSE
149:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr);

151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: !  Create nonlinear solver context
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

155:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

157:       if (pc .eqv. PETSC_TRUE) then
158:         call SNESSetType(snes,SNESNEWTONTR,ierr)
159:         call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)
160:       endif

162: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: !  Create vector data structures; set function evaluation routine
164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

166:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
167:       call VecSetSizes(x,PETSC_DECIDE,N,ierr)
168:       call VecSetFromOptions(x,ierr)
169:       call VecDuplicate(x,r,ierr)

171: !  Set function evaluation routine and vector.  Whenever the nonlinear
172: !  solver needs to evaluate the nonlinear function, it will call this
173: !  routine.
174: !   - Note that the final routine argument is the user-defined
175: !     context that provides application-specific data for the
176: !     function evaluation routine.

178:       call SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr)

180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: !  Create matrix data structure; set Jacobian evaluation routine
182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

188:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)
189:       if (.not. matrix_free) then
190:         call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)
191:       endif

193: !
194: !     This option will cause the Jacobian to be computed via finite differences
195: !    efficiently using a coloring of the columns of the matrix.
196: !
197:       fd_coloring = .false.
198:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr)
199:       if (fd_coloring) then

201: !
202: !      This initializes the nonzero structure of the Jacobian. This is artificial
203: !      because clearly if we had a routine to compute the Jacobian we won't need
204: !      to use finite differences.
205: !
206:         call FormJacobian(snes,x,J,J,0,ierr)
207: !
208: !       Color the matrix, i.e. determine groups of columns that share no common
209: !      rows. These columns in the Jacobian can all be computed simultaneously.
210: !
211:         call MatColoringCreate(J,mc,ierr)
212:         call MatColoringSetType(mc,MATCOLORINGNATURAL,ierr)
213:         call MatColoringSetFromOptions(mc,ierr)
214:         call MatColoringApply(mc,iscoloring,ierr)
215:         call MatColoringDestroy(mc,ierr)
216: !
217: !       Create the data structure that SNESComputeJacobianDefaultColor() uses
218: !       to compute the actual Jacobians via finite differences.
219: !
220:         call MatFDColoringCreate(J,iscoloring,fdcoloring,ierr)
221:         call MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr)
222:         call MatFDColoringSetFromOptions(fdcoloring,ierr)
223:         call MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr)
224: !
225: !        Tell SNES to use the routine SNESComputeJacobianDefaultColor()
226: !      to compute Jacobians.
227: !
228:         call SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr)
229:         call ISColoringDestroy(iscoloring,ierr)

231:       else if (.not. matrix_free) then

233: !  Set Jacobian matrix data structure and default Jacobian evaluation
234: !  routine.  Whenever the nonlinear solver needs to compute the
235: !  Jacobian matrix, it will call this routine.
236: !   - Note that the final routine argument is the user-defined
237: !     context that provides application-specific data for the
238: !     Jacobian evaluation routine.
239: !   - The user can override with:
240: !      -snes_fd : default finite differencing approximation of Jacobian
241: !      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
242: !                 (unless user explicitly sets preconditioner)
243: !      -snes_mf_operator : form preconditioning matrix as set by the user,
244: !                          but use matrix-free approx for Jacobian-vector
245: !                          products within Newton-Krylov method
246: !
247:         call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
248:       endif

250: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251: !  Customize nonlinear solver; set runtime options
252: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

256:       call SNESSetFromOptions(snes,ierr)

258: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259: !  Evaluate initial guess; then solve nonlinear system.
260: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

267:       call FormInitialGuess(x,ierr)
268:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
269:       call SNESGetIterationNumber(snes,its,ierr);
270:       if (rank .eq. 0) then
271:          write(6,100) its
272:       endif
273:   100 format('Number of SNES iterations = ',i1)

275: !  PetscDraw contour plot of solution

277:       call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
278:       call PetscDrawSetFromOptions(draw,ierr)

280:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
281:       call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)
282:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

284: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
285: !  Free work space.  All PETSc objects should be destroyed when they
286: !  are no longer needed.
287: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

289:       if (.not. matrix_free) call MatDestroy(J,ierr)
290:       if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)

292:       call VecDestroy(x,ierr)
293:       call VecDestroy(r,ierr)
294:       call SNESDestroy(snes,ierr)
295:       call PetscDrawDestroy(draw,ierr)
296:       call PetscFinalize(ierr)
297:       end

299: ! ---------------------------------------------------------------------
300: !
301: !  FormInitialGuess - Forms initial approximation.
302: !
303: !  Input Parameter:
304: !  X - vector
305: !
306: !  Output Parameters:
307: !  X - vector
308: !  ierr - error code
309: !
310: !  Notes:
311: !  This routine serves as a wrapper for the lower-level routine
312: !  "ApplicationInitialGuess", where the actual computations are
313: !  done using the standard Fortran style of treating the local
314: !  vector data as a multidimensional array over the local mesh.
315: !  This routine merely accesses the local vector data via
316: !  VecGetArray() and VecRestoreArray().
317: !
318:       subroutine FormInitialGuess(X,ierr)
319:       use petscsnes
320:       implicit none

322: !  Input/output variables:
323:       Vec           X
324:       PetscErrorCode    ierr

326: !  Declarations for use with local arrays:
327:       PetscScalar   lx_v(0:1)
328:       PetscOffset   lx_i

330:       0

332: !  Get a pointer to vector data.
333: !    - For default PETSc vectors, VecGetArray() returns a pointer to
334: !      the data array.  Otherwise, the routine is implementation dependent.
335: !    - You MUST call VecRestoreArray() when you no longer need access to
336: !      the array.
337: !    - Note that the Fortran interface to VecGetArray() differs from the
338: !      C version.  See the users manual for details.

340:       call VecGetArray(X,lx_v,lx_i,ierr)

342: !  Compute initial guess

344:       call ApplicationInitialGuess(lx_v(lx_i),ierr)

346: !  Restore vector

348:       call VecRestoreArray(X,lx_v,lx_i,ierr)

350:       return
351:       end

353: ! ---------------------------------------------------------------------
354: !
355: !  ApplicationInitialGuess - Computes initial approximation, called by
356: !  the higher level routine FormInitialGuess().
357: !
358: !  Input Parameter:
359: !  x - local vector data
360: !
361: !  Output Parameters:
362: !  f - local vector data, f(x)
363: !  ierr - error code
364: !
365: !  Notes:
366: !  This routine uses standard Fortran-style computations over a 2-dim array.
367: !
368:       subroutine ApplicationInitialGuess(x,ierr)
369:       use petscksp
370:       implicit none

372: !  Common blocks:
373:       PetscReal   lambda
374:       PetscInt     mx,my
375:       PetscBool         fd_coloring
376:       common      /params/ lambda,mx,my,fd_coloring

378: !  Input/output variables:
379:       PetscScalar x(mx,my)
380:       PetscErrorCode     ierr

382: !  Local variables:
383:       PetscInt     i,j
384:       PetscReal temp1,temp,hx,hy,one

386: !  Set parameters

388:       0
389:       one    = 1.0
390:       hx     = one/(mx-1)
391:       hy     = one/(my-1)
392:       temp1  = lambda/(lambda + one)

394:       do 20 j=1,my
395:          temp = min(j-1,my-j)*hy
396:          do 10 i=1,mx
397:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
398:               x(i,j) = 0.0
399:             else
400:               x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
401:             endif
402:  10      continue
403:  20   continue

405:       return
406:       end

408: ! ---------------------------------------------------------------------
409: !
410: !  FormFunction - Evaluates nonlinear function, F(x).
411: !
412: !  Input Parameters:
413: !  snes  - the SNES context
414: !  X     - input vector
415: !  dummy - optional user-defined context, as set by SNESSetFunction()
416: !          (not used here)
417: !
418: !  Output Parameter:
419: !  F     - vector with newly computed function
420: !
421: !  Notes:
422: !  This routine serves as a wrapper for the lower-level routine
423: !  "ApplicationFunction", where the actual computations are
424: !  done using the standard Fortran style of treating the local
425: !  vector data as a multidimensional array over the local mesh.
426: !  This routine merely accesses the local vector data via
427: !  VecGetArray() and VecRestoreArray().
428: !
429:       subroutine FormFunction(snes,X,F,fdcoloring,ierr)
430:       use petscsnes
431:       implicit none

433: !  Input/output variables:
434:       SNES              snes
435:       Vec               X,F
436:       PetscErrorCode          ierr
437:       MatFDColoring fdcoloring

439: !  Common blocks:
440:       PetscReal         lambda
441:       PetscInt          mx,my
442:       PetscBool         fd_coloring
443:       common            /params/ lambda,mx,my,fd_coloring

445: !  Declarations for use with local arrays:
446:       PetscScalar       lx_v(0:1),lf_v(0:1)
447:       PetscOffset       lx_i,lf_i

449:       PetscInt, pointer :: indices(:)

451: !  Get pointers to vector data.
452: !    - For default PETSc vectors, VecGetArray() returns a pointer to
453: !      the data array.  Otherwise, the routine is implementation dependent.
454: !    - You MUST call VecRestoreArray() when you no longer need access to
455: !      the array.
456: !    - Note that the Fortran interface to VecGetArray() differs from the
457: !      C version.  See the Fortran chapter of the users manual for details.

459:       call VecGetArrayRead(X,lx_v,lx_i,ierr)
460:       call VecGetArray(F,lf_v,lf_i,ierr)

462: !  Compute function

464:       call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)

466: !  Restore vectors

468:       call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
469:       call VecRestoreArray(F,lf_v,lf_i,ierr)

471:       call PetscLogFlops(11.0d0*mx*my,ierr)
472: !
473: !     fdcoloring is in the common block and used here ONLY to test the
474: !     calls to MatFDColoringGetPerturbedColumnsF90() and  MatFDColoringRestorePerturbedColumnsF90()
475: !
476:       if (fd_coloring) then
477:          call MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)
478:          print*,'Indices from GetPerturbedColumnsF90'
479:          write(*,1000) indices
480:  1000    format(50i4)
481:          call MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)
482:       endif
483:       return
484:       end

486: ! ---------------------------------------------------------------------
487: !
488: !  ApplicationFunction - Computes nonlinear function, called by
489: !  the higher level routine FormFunction().
490: !
491: !  Input Parameter:
492: !  x    - local vector data
493: !
494: !  Output Parameters:
495: !  f    - local vector data, f(x)
496: !  ierr - error code
497: !
498: !  Notes:
499: !  This routine uses standard Fortran-style computations over a 2-dim array.
500: !
501:       subroutine ApplicationFunction(x,f,ierr)
502:       use petscsnes
503:       implicit none

505: !  Common blocks:
506:       PetscReal      lambda
507:       PetscInt        mx,my
508:       PetscBool         fd_coloring
509:       common         /params/ lambda,mx,my,fd_coloring

511: !  Input/output variables:
512:       PetscScalar    x(mx,my),f(mx,my)
513:       PetscErrorCode       ierr

515: !  Local variables:
516:       PetscScalar    two,one,hx,hy
517:       PetscScalar    hxdhy,hydhx,sc
518:       PetscScalar    u,uxx,uyy
519:       PetscInt        i,j

521:       0
522:       one    = 1.0
523:       two    = 2.0
524:       hx     = one/(mx-1)
525:       hy     = one/(my-1)
526:       sc     = hx*hy*lambda
527:       hxdhy  = hx/hy
528:       hydhx  = hy/hx

530: !  Compute function

532:       do 20 j=1,my
533:          do 10 i=1,mx
534:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
535:                f(i,j) = x(i,j)
536:             else
537:                u = x(i,j)
538:                uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
539:                uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
540:                f(i,j) = uxx + uyy - sc*exp(u)
541:             endif
542:  10      continue
543:  20   continue

545:       return
546:       end

548: ! ---------------------------------------------------------------------
549: !
550: !  FormJacobian - Evaluates Jacobian matrix.
551: !
552: !  Input Parameters:
553: !  snes    - the SNES context
554: !  x       - input vector
555: !  dummy   - optional user-defined context, as set by SNESSetJacobian()
556: !            (not used here)
557: !
558: !  Output Parameters:
559: !  jac      - Jacobian matrix
560: !  jac_prec - optionally different preconditioning matrix (not used here)
561: !  flag     - flag indicating matrix structure
562: !
563: !  Notes:
564: !  This routine serves as a wrapper for the lower-level routine
565: !  "ApplicationJacobian", where the actual computations are
566: !  done using the standard Fortran style of treating the local
567: !  vector data as a multidimensional array over the local mesh.
568: !  This routine merely accesses the local vector data via
569: !  VecGetArray() and VecRestoreArray().
570: !
571:       subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
572:       use petscsnes
573:       implicit none

575: !  Input/output variables:
576:       SNES          snes
577:       Vec           X
578:       Mat           jac,jac_prec
579:       PetscErrorCode      ierr
580:       integer dummy

582: !  Common blocks:
583:       PetscReal     lambda
584:       PetscInt       mx,my
585:       PetscBool         fd_coloring
586:       common        /params/ lambda,mx,my,fd_coloring

588: !  Declarations for use with local array:
589:       PetscScalar   lx_v(0:1)
590:       PetscOffset   lx_i

592: !  Get a pointer to vector data

594:       call VecGetArrayRead(X,lx_v,lx_i,ierr)

596: !  Compute Jacobian entries

598:       call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)

600: !  Restore vector

602:       call VecRestoreArrayRead(X,lx_v,lx_i,ierr)

604: !  Assemble matrix

606:       call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
607:       call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)

609:       return
610:       end

612: ! ---------------------------------------------------------------------
613: !
614: !  ApplicationJacobian - Computes Jacobian matrix, called by
615: !  the higher level routine FormJacobian().
616: !
617: !  Input Parameters:
618: !  x        - local vector data
619: !
620: !  Output Parameters:
621: !  jac      - Jacobian matrix
622: !  jac_prec - optionally different preconditioning matrix (not used here)
623: !  ierr     - error code
624: !
625: !  Notes:
626: !  This routine uses standard Fortran-style computations over a 2-dim array.
627: !
628:       subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
629:       use petscsnes
630:       implicit none

632: !  Common blocks:
633:       PetscReal    lambda
634:       PetscInt      mx,my
635:       PetscBool         fd_coloring
636:       common       /params/ lambda,mx,my,fd_coloring

638: !  Input/output variables:
639:       PetscScalar  x(mx,my)
640:       Mat          jac,jac_prec
641:       PetscErrorCode      ierr

643: !  Local variables:
644:       PetscInt      i,j,row(1),col(5),i1,i5
645:       PetscScalar  two,one, hx,hy
646:       PetscScalar  hxdhy,hydhx,sc,v(5)

648: !  Set parameters

650:       i1     = 1
651:       i5     = 5
652:       one    = 1.0
653:       two    = 2.0
654:       hx     = one/(mx-1)
655:       hy     = one/(my-1)
656:       sc     = hx*hy
657:       hxdhy  = hx/hy
658:       hydhx  = hy/hx

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

665:       do 20 j=1,my
666:          row(1) = (j-1)*mx - 1
667:          do 10 i=1,mx
668:             row(1) = row(1) + 1
669: !           boundary points
670:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
671:                call MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr)
672: !           interior grid points
673:             else
674:                v(1) = -hxdhy
675:                v(2) = -hydhx
676:                v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
677:                v(4) = -hydhx
678:                v(5) = -hxdhy
679:                col(1) = row(1) - mx
680:                col(2) = row(1) - 1
681:                col(3) = row(1)
682:                col(4) = row(1) + 1
683:                col(5) = row(1) + mx
684:                call MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr)
685:             endif
686:  10      continue
687:  20   continue

689:       return
690:       end

692: !
693: !/*TEST
694: !
695: !   build:
696: !      requires: !single
697: !
698: !   test:
699: !      args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
700: !
701: !   test:
702: !      suffix: 2
703: !      args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
704: !
705: !   test:
706: !      suffix: 3
707: !      args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
708: !      filter: sort -b
709: !      filter_output: sort -b
710: !
711: !   test:
712: !     suffix: 4
713: !     args: -pc -par 6.807 -nox
714: !
715: !TEST*/