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: subroutine postcheck(snes,x,y,w,changed_y,changed_w,ctx,ierr)
31: #include <petsc/finclude/petscsnes.h>
32: use petscsnes
33: implicit none
34: SNES snes
35: PetscReal norm
36: Vec tmp,x,y,w
37: PetscBool changed_w,changed_y
38: PetscErrorCode ierr
39: PetscInt ctx
40: PetscScalar mone
42: PetscCallA(VecDuplicate(x,tmp,ierr))
43: mone = -1.0
44: PetscCallA(VecWAXPY(tmp,mone,x,w,ierr))
45: PetscCallA(VecNorm(tmp,NORM_2,norm,ierr))
46: PetscCallA(VecDestroy(tmp,ierr))
47: print*, 'Norm of search step ',norm
48: changed_y = PETSC_FALSE
49: changed_w = PETSC_FALSE
50: return
51: end
53: program main
54: #include <petsc/finclude/petscdraw.h>
55: use petscsnes
56: implicit none
57: interface SNESSetJacobian
58: subroutine SNESSetJacobian1(a,b,c,d,e,z)
59: use petscsnes
60: SNES a
61: Mat b
62: Mat c
63: external d
64: MatFDColoring e
65: PetscErrorCode z
66: end subroutine
67: subroutine SNESSetJacobian2(a,b,c,d,e,z)
68: use petscsnes
69: SNES a
70: Mat b
71: Mat c
72: external d
73: integer e
74: PetscErrorCode z
75: end subroutine
76: end interface
77: !
78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: ! Variable declarations
80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: !
82: ! Variables:
83: ! snes - nonlinear solver
84: ! x, r - solution, residual vectors
85: ! J - Jacobian matrix
86: ! its - iterations for convergence
87: ! matrix_free - flag - 1 indicates matrix-free version
88: ! lambda - nonlinearity parameter
89: ! draw - drawing context
90: !
91: SNES snes
92: MatColoring mc
93: Vec x,r
94: PetscDraw draw
95: Mat J
96: PetscBool matrix_free,flg,fd_coloring
97: PetscErrorCode ierr
98: PetscInt its,N, mx,my,i5
99: PetscMPIInt size,rank
100: PetscReal lambda_max,lambda_min,lambda
101: MatFDColoring fdcoloring
102: ISColoring iscoloring
103: PetscBool pc
104: external postcheck
106: PetscScalar,pointer :: lx_v(:)
108: ! Store parameters in common block
110: common /params/ lambda,mx,my,fd_coloring
112: ! Note: Any user-defined Fortran routines (such as FormJacobian)
113: ! MUST be declared as external.
115: external FormFunction,FormInitialGuess,FormJacobian
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: ! Initialize program
119: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: PetscCallA(PetscInitialize(ierr))
122: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
123: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
125: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
127: ! Initialize problem parameters
128: i5 = 5
129: lambda_max = 6.81
130: lambda_min = 0.0
131: lambda = 6.0
132: mx = 4
133: my = 4
134: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
135: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
136: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr))
137: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range '); endif
138: N = mx*my
139: pc = PETSC_FALSE
140: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr))
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: ! Create nonlinear solver context
144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))
148: if (pc .eqv. PETSC_TRUE) then
149: PetscCallA(SNESSetType(snes,SNESNEWTONTR,ierr))
150: PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr))
151: endif
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: ! Create vector data structures; set function evaluation routine
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
158: PetscCallA(VecSetSizes(x,PETSC_DECIDE,N,ierr))
159: PetscCallA(VecSetFromOptions(x,ierr))
160: PetscCallA(VecDuplicate(x,r,ierr))
162: ! Set function evaluation routine and vector. Whenever the nonlinear
163: ! solver needs to evaluate the nonlinear function, it will call this
164: ! routine.
165: ! - Note that the final routine argument is the user-defined
166: ! context that provides application-specific data for the
167: ! function evaluation routine.
169: PetscCallA(SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr))
171: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: ! Create matrix data structure; set Jacobian evaluation routine
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: ! Create matrix. Here we only approximately preallocate storage space
176: ! for the Jacobian. See the users manual for a discussion of better
177: ! techniques for preallocating matrix memory.
179: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr))
180: if (.not. matrix_free) then
181: PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr))
182: endif
184: !
185: ! This option will cause the Jacobian to be computed via finite differences
186: ! efficiently using a coloring of the columns of the matrix.
187: !
188: fd_coloring = .false.
189: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr))
190: if (fd_coloring) then
192: !
193: ! This initializes the nonzero structure of the Jacobian. This is artificial
194: ! because clearly if we had a routine to compute the Jacobian we won't need
195: ! to use finite differences.
196: !
197: PetscCallA(FormJacobian(snes,x,J,J,0,ierr))
198: !
199: ! Color the matrix, i.e. determine groups of columns that share no common
200: ! rows. These columns in the Jacobian can all be computed simultaneously.
201: !
202: PetscCallA(MatColoringCreate(J,mc,ierr))
203: PetscCallA(MatColoringSetType(mc,MATCOLORINGNATURAL,ierr))
204: PetscCallA(MatColoringSetFromOptions(mc,ierr))
205: PetscCallA(MatColoringApply(mc,iscoloring,ierr))
206: PetscCallA(MatColoringDestroy(mc,ierr))
207: !
208: ! Create the data structure that SNESComputeJacobianDefaultColor() uses
209: ! to compute the actual Jacobians via finite differences.
210: !
211: PetscCallA(MatFDColoringCreate(J,iscoloring,fdcoloring,ierr))
212: PetscCallA(MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr))
213: PetscCallA(MatFDColoringSetFromOptions(fdcoloring,ierr))
214: PetscCallA(MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr))
215: !
216: ! Tell SNES to use the routine SNESComputeJacobianDefaultColor()
217: ! to compute Jacobians.
218: !
219: PetscCallA(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr))
220: PetscCallA(ISColoringDestroy(iscoloring,ierr))
222: else if (.not. matrix_free) then
224: ! Set Jacobian matrix data structure and default Jacobian evaluation
225: ! routine. Whenever the nonlinear solver needs to compute the
226: ! Jacobian matrix, it will call this routine.
227: ! - Note that the final routine argument is the user-defined
228: ! context that provides application-specific data for the
229: ! Jacobian evaluation routine.
230: ! - The user can override with:
231: ! -snes_fd : default finite differencing approximation of Jacobian
232: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
233: ! (unless user explicitly sets preconditioner)
234: ! -snes_mf_operator : form preconditioning matrix as set by the user,
235: ! but use matrix-free approx for Jacobian-vector
236: ! products within Newton-Krylov method
237: !
238: PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
239: endif
241: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
242: ! Customize nonlinear solver; set runtime options
243: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
245: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
247: PetscCallA(SNESSetFromOptions(snes,ierr))
249: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250: ! Evaluate initial guess; then solve nonlinear system.
251: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253: ! Note: The user should initialize the vector, x, with the initial guess
254: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
255: ! to employ an initial guess of zero, the user should explicitly set
256: ! this vector to zero by calling VecSet().
258: PetscCallA(FormInitialGuess(x,ierr))
259: PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
260: PetscCallA(SNESGetIterationNumber(snes,its,ierr))
261: if (rank .eq. 0) then
262: write(6,100) its
263: endif
264: 100 format('Number of SNES iterations = ',i1)
266: ! PetscDraw contour plot of solution
268: PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr))
269: PetscCallA(PetscDrawSetFromOptions(draw,ierr))
271: PetscCallA(VecGetArrayReadF90(x,lx_v,ierr))
272: PetscCallA(PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v,ierr))
273: PetscCallA(VecRestoreArrayReadF90(x,lx_v,ierr))
275: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276: ! Free work space. All PETSc objects should be destroyed when they
277: ! are no longer needed.
278: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280: if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr))
281: if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring,ierr))
283: PetscCallA(VecDestroy(x,ierr))
284: PetscCallA(VecDestroy(r,ierr))
285: PetscCallA(SNESDestroy(snes,ierr))
286: PetscCallA(PetscDrawDestroy(draw,ierr))
287: PetscCallA(PetscFinalize(ierr))
288: end
290: ! ---------------------------------------------------------------------
291: !
292: ! FormInitialGuess - Forms initial approximation.
293: !
294: ! Input Parameter:
295: ! X - vector
296: !
297: ! Output Parameters:
298: ! X - vector
299: ! ierr - error code
300: !
301: ! Notes:
302: ! This routine serves as a wrapper for the lower-level routine
303: ! "ApplicationInitialGuess", where the actual computations are
304: ! done using the standard Fortran style of treating the local
305: ! vector data as a multidimensional array over the local mesh.
306: ! This routine merely accesses the local vector data via
307: ! VecGetArrayF90() and VecRestoreArrayF90().
308: !
309: subroutine FormInitialGuess(X,ierr)
310: use petscsnes
311: implicit none
313: ! Input/output variables:
314: Vec X
315: PetscErrorCode ierr
317: ! Declarations for use with local arrays:
318: PetscScalar,pointer :: lx_v(:)
320: ierr = 0
322: ! Get a pointer to vector data.
323: ! - VecGetArrayF90() returns a pointer to the data array.
324: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
325: ! the array.
327: PetscCallA(VecGetArrayF90(X,lx_v,ierr))
329: ! Compute initial guess
331: PetscCallA(ApplicationInitialGuess(lx_v,ierr))
333: ! Restore vector
335: PetscCallA(VecRestoreArrayF90(X,lx_v,ierr))
337: return
338: end
340: ! ApplicationInitialGuess - Computes initial approximation, called by
341: ! the higher level routine FormInitialGuess().
342: !
343: ! Input Parameter:
344: ! x - local vector data
345: !
346: ! Output Parameters:
347: ! f - local vector data, f(x)
348: ! ierr - error code
349: !
350: ! Notes:
351: ! This routine uses standard Fortran-style computations over a 2-dim array.
352: !
353: subroutine ApplicationInitialGuess(x,ierr)
354: use petscksp
355: implicit none
357: ! Common blocks:
358: PetscReal lambda
359: PetscInt mx,my
360: PetscBool fd_coloring
361: common /params/ lambda,mx,my,fd_coloring
363: ! Input/output variables:
364: PetscScalar x(mx,my)
365: PetscErrorCode ierr
367: ! Local variables:
368: PetscInt i,j
369: PetscReal temp1,temp,hx,hy,one
371: ! Set parameters
373: ierr = 0
374: one = 1.0
375: hx = one/(mx-1)
376: hy = one/(my-1)
377: temp1 = lambda/(lambda + one)
379: do 20 j=1,my
380: temp = min(j-1,my-j)*hy
381: do 10 i=1,mx
382: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
383: x(i,j) = 0.0
384: else
385: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
386: endif
387: 10 continue
388: 20 continue
390: return
391: end
393: ! ---------------------------------------------------------------------
394: !
395: ! FormFunction - Evaluates nonlinear function, F(x).
396: !
397: ! Input Parameters:
398: ! snes - the SNES context
399: ! X - input vector
400: ! dummy - optional user-defined context, as set by SNESSetFunction()
401: ! (not used here)
402: !
403: ! Output Parameter:
404: ! F - vector with newly computed function
405: !
406: ! Notes:
407: ! This routine serves as a wrapper for the lower-level routine
408: ! "ApplicationFunction", where the actual computations are
409: ! done using the standard Fortran style of treating the local
410: ! vector data as a multidimensional array over the local mesh.
411: ! This routine merely accesses the local vector data via
412: ! VecGetArrayF90() and VecRestoreArrayF90().
413: !
414: subroutine FormFunction(snes,X,F,fdcoloring,ierr)
415: use petscsnes
416: implicit none
418: ! Input/output variables:
419: SNES snes
420: Vec X,F
421: PetscErrorCode ierr
422: MatFDColoring fdcoloring
424: ! Common blocks:
425: PetscReal lambda
426: PetscInt mx,my
427: PetscBool fd_coloring
428: common /params/ lambda,mx,my,fd_coloring
430: ! Declarations for use with local arrays:
431: PetscScalar,pointer :: lx_v(:), lf_v(:)
432: PetscInt, pointer :: indices(:)
434: ! Get pointers to vector data.
435: ! - VecGetArrayF90() returns a pointer to the data array.
436: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
437: ! the array.
439: PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))
440: PetscCallA(VecGetArrayF90(F,lf_v,ierr))
442: ! Compute function
444: PetscCallA(ApplicationFunction(lx_v,lf_v,ierr))
446: ! Restore vectors
448: PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))
449: PetscCallA(VecRestoreArrayF90(F,lf_v,ierr))
451: PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr))
452: !
453: ! fdcoloring is in the common block and used here ONLY to test the
454: ! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90()
455: !
456: if (fd_coloring) then
457: PetscCallA(MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr))
458: print*,'Indices from GetPerturbedColumnsF90'
459: write(*,1000) indices
460: 1000 format(50i4)
461: PetscCallA(MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr))
462: endif
463: return
464: end
466: ! ---------------------------------------------------------------------
467: !
468: ! ApplicationFunction - Computes nonlinear function, called by
469: ! the higher level routine FormFunction().
470: !
471: ! Input Parameter:
472: ! x - local vector data
473: !
474: ! Output Parameters:
475: ! f - local vector data, f(x)
476: ! ierr - error code
477: !
478: ! Notes:
479: ! This routine uses standard Fortran-style computations over a 2-dim array.
480: !
481: subroutine ApplicationFunction(x,f,ierr)
482: use petscsnes
483: implicit none
485: ! Common blocks:
486: PetscReal lambda
487: PetscInt mx,my
488: PetscBool fd_coloring
489: common /params/ lambda,mx,my,fd_coloring
491: ! Input/output variables:
492: PetscScalar x(mx,my),f(mx,my)
493: PetscErrorCode ierr
495: ! Local variables:
496: PetscScalar two,one,hx,hy
497: PetscScalar hxdhy,hydhx,sc
498: PetscScalar u,uxx,uyy
499: PetscInt i,j
501: ierr = 0
502: one = 1.0
503: two = 2.0
504: hx = one/(mx-1)
505: hy = one/(my-1)
506: sc = hx*hy*lambda
507: hxdhy = hx/hy
508: hydhx = hy/hx
510: ! Compute function
512: do 20 j=1,my
513: do 10 i=1,mx
514: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
515: f(i,j) = x(i,j)
516: else
517: u = x(i,j)
518: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
519: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
520: f(i,j) = uxx + uyy - sc*exp(u)
521: endif
522: 10 continue
523: 20 continue
525: return
526: end
528: ! ---------------------------------------------------------------------
529: !
530: ! FormJacobian - Evaluates Jacobian matrix.
531: !
532: ! Input Parameters:
533: ! snes - the SNES context
534: ! x - input vector
535: ! dummy - optional user-defined context, as set by SNESSetJacobian()
536: ! (not used here)
537: !
538: ! Output Parameters:
539: ! jac - Jacobian matrix
540: ! jac_prec - optionally different preconditioning matrix (not used here)
541: ! flag - flag indicating matrix structure
542: !
543: ! Notes:
544: ! This routine serves as a wrapper for the lower-level routine
545: ! "ApplicationJacobian", where the actual computations are
546: ! done using the standard Fortran style of treating the local
547: ! vector data as a multidimensional array over the local mesh.
548: ! This routine merely accesses the local vector data via
549: ! VecGetArrayF90() and VecRestoreArrayF90().
550: !
551: subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
552: use petscsnes
553: implicit none
555: ! Input/output variables:
556: SNES snes
557: Vec X
558: Mat jac,jac_prec
559: PetscErrorCode ierr
560: integer dummy
562: ! Common blocks:
563: PetscReal lambda
564: PetscInt mx,my
565: PetscBool fd_coloring
566: common /params/ lambda,mx,my,fd_coloring
568: ! Declarations for use with local array:
569: PetscScalar,pointer :: lx_v(:)
571: ! Get a pointer to vector data
573: PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))
575: ! Compute Jacobian entries
577: PetscCallA(ApplicationJacobian(lx_v,jac,jac_prec,ierr))
579: ! Restore vector
581: PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))
583: ! Assemble matrix
585: PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
586: PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
588: return
589: end
591: ! ---------------------------------------------------------------------
592: !
593: ! ApplicationJacobian - Computes Jacobian matrix, called by
594: ! the higher level routine FormJacobian().
595: !
596: ! Input Parameters:
597: ! x - local vector data
598: !
599: ! Output Parameters:
600: ! jac - Jacobian matrix
601: ! jac_prec - optionally different preconditioning matrix (not used here)
602: ! ierr - error code
603: !
604: ! Notes:
605: ! This routine uses standard Fortran-style computations over a 2-dim array.
606: !
607: subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
608: use petscsnes
609: implicit none
611: ! Common blocks:
612: PetscReal lambda
613: PetscInt mx,my
614: PetscBool fd_coloring
615: common /params/ lambda,mx,my,fd_coloring
617: ! Input/output variables:
618: PetscScalar x(mx,my)
619: Mat jac,jac_prec
620: PetscErrorCode ierr
622: ! Local variables:
623: PetscInt i,j,row(1),col(5),i1,i5
624: PetscScalar two,one, hx,hy
625: PetscScalar hxdhy,hydhx,sc,v(5)
627: ! Set parameters
629: i1 = 1
630: i5 = 5
631: one = 1.0
632: two = 2.0
633: hx = one/(mx-1)
634: hy = one/(my-1)
635: sc = hx*hy
636: hxdhy = hx/hy
637: hydhx = hy/hx
639: ! Compute entries of the Jacobian matrix
640: ! - Here, we set all entries for a particular row at once.
641: ! - Note that MatSetValues() uses 0-based row and column numbers
642: ! in Fortran as well as in C.
644: do 20 j=1,my
645: row(1) = (j-1)*mx - 1
646: do 10 i=1,mx
647: row(1) = row(1) + 1
648: ! boundary points
649: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
650: PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr))
651: ! interior grid points
652: else
653: v(1) = -hxdhy
654: v(2) = -hydhx
655: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
656: v(4) = -hydhx
657: v(5) = -hxdhy
658: col(1) = row(1) - mx
659: col(2) = row(1) - 1
660: col(3) = row(1)
661: col(4) = row(1) + 1
662: col(5) = row(1) + mx
663: PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr))
664: endif
665: 10 continue
666: 20 continue
668: return
669: end
671: !
672: !/*TEST
673: !
674: ! build:
675: ! requires: !single
676: !
677: ! test:
678: ! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
679: !
680: ! test:
681: ! suffix: 2
682: ! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
683: !
684: ! test:
685: ! suffix: 3
686: ! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
687: ! filter: sort -b
688: ! filter_output: sort -b
689: !
690: ! test:
691: ! suffix: 4
692: ! args: -pc -par 6.807 -nox
693: !
694: !TEST*/